Upper Colorado River basin streamflow has declined by roughly 20% over the last century of the instrumental period, based on estimates of naturalized flow above Lees Ferry. Here we assess factors causing the decline and evaluate the premise that rising surface temperatures have been mostly responsible. We use an event attribution framework involving parallel sets of global model experiments with and without climate change drivers. We demonstrate that climate change forcing has acted to reduce Upper Colorado River basin streamflow during this period by about 10% (with uncertainty range of 6%–14% reductions). The magnitude of the observed flow decline is found to be inconsistent with natural variability alone, and approximately one-half of the observed flow decline is judged to have resulted from long-term climate change. Each of three different global models used herein indicates that climate change forcing during the last century has acted to increase surface temperature (~+1.2°C) and decrease precipitation (~−3%). Using large ensemble methods, we diagnose the separate effects of temperature and precipitation changes on Upper Colorado River streamflow. Precipitation change is found to be the most consequential factor owing to its amplified impact on flow resulting from precipitation elasticity (percent change in streamflow per percent change in precipitation) of ~2. We confirm that warming has also driven streamflow declines, as inferred from empirical studies, although operating as a secondary factor. Our finding of a modest −2.5% °C−1 temperature sensitivity, on the basis of our best model-derived estimate, indicates that only about one-third of the attributable climate change signal in Colorado River decline resulted from warming, whereas about two-thirds resulted from precipitation decline.
Colorado River streamflow has fallen substantially during the period of the instrumental record. A linear trend of Upper Colorado River Commission naturalized Lees Ferry flow over 1896–2018 indicates a 20% decline (Fig. 1). For 1906–2018, a span over which the U.S. Bureau of Reclamation has also produced a naturalized flow record, a 24% decline is indicated by both datasets. The Lees Ferry stream gauge measures surface water flow accumulated from tributaries draining the Upper Colorado River basin (UCRB), and this upper reach contains the principal headwaters that account for almost 90% of total Colorado River basin streamflow (Jacobs 2011). Its decrease, coupled with recent demand for water eclipsing supply (Rajagopalan et al. 2009), is thus of concern to large parts of the southwestern United States and Mexico sharing water under the Law of the River for which the Colorado River Compact is cornerstone (https://www.usbr.gov/lc/region/g1000/lawofrvr.html).
One hypothesis for secular decline in the instrumental period1—initially proposed 40 years ago (Stockton and Boggess 1979; Revelle and Waggoner 1983) and recently resurrected (Hoerling and Eischeid 2007; Woodhouse et al. 2016; Udall and Overpeck 2017; McCabe et al. 2017)—holds that Colorado River streamflow is highly sensitive to temperature such that warming appreciably reduces flow volumes. A warming over the last century (Fig. 1) in the UCRB (Hoerling et al. 2013; Lukas et al. 2014) has itself been linked to anthropogenic climate change, inducing aridification across the greater U.S. Southwest (Seager et al. 2007; Cook et al. 2015). A hypothesis that argues for strong temperature effects on historical Colorado River flow is thus a portent for future severe declines owing to high confidence that warming will continue to occur (IPCC 2013).
Guidance on future Colorado River flows thus rests heavily on explaining its decline in the instrumental period, and especially on determining its sensitivity to temperature. However, the magnitude of that sensitivity remains highly uncertain despite a relatively long history of studies on the problem mentioned above. Wide ranging estimates exist that speak to complications in deciphering empirical results and to difficulties in interpreting model-based studies. For instance, different land surface model (LSM) simulations indicate Colorado River flow might decline as little as 3% or perhaps as much as 8% per degree Celsius of warming (Christensen et al. 2004; Christensen and Lettenmaier 2007; Vano et al. 2012, 2014). Empirical studies suggest sensitivities from −10% to −15% °C−1 (Revelle and Waggoner 1983; Nowak et al. 2012; McCabe et al. 2017), although sparse observations and sampling variability lead to large error bars in attributing effects of meteorological drivers (temperature and precipitation) on streamflow variability (Milly and Dunne 2002). In a different approach to this problem, an extension of the theory encapsulated by Budyko (1974) on climatological water and energy balances indicates low sensitivity of streamflow to temperature, almost never exceeding about −4% °C−1 (Milly et al. 2018). Suffice it to say that expectations for severe late twenty-first-century declines in Colorado River flow would be warranted if sensitivity to temperature reaches the high end of these estimates.
Here, we apply new methods to test hypotheses of causality and especially the proposition for a large temperature role in Colorado River–climate linkages. We employ an event attribution framework that has been previously used to quantify climate change effects on weather extremes (National Academies of Sciences, Engineering, and Medicine 2016). As described in section 2, this involves a dynamical modeling approach whereby large ensembles of simulations from multiple atmosphere/land models are subjected to various drivers selected to distinguish effects of climate change from natural variability. The approach is in analogy with epidemiology wherein determinants of disease in human populations are of interest, and possible cures are evaluated by comparing two population samples, with one subjected to a potential curing drug and the other given a placebo. In our hydroclimate study, the placebo population consists of experiments without global warming influence, referred to as counterfactual experiments, whose statistics are compared with a second set of runs incorporating all current forcings (factual experiments). Our model experiments are conducted globally, recognizing that UCRB hydrology responds to worldwide drivers such as naturally occurring Pacific and Atlantic Ocean variability (McCabe et al. 2004; Nowak et al. 2012) and to global climate change (Cayan et al. 2010; Seager and Vecchi 2010).
We address various limitations in prior modeling investigations of Colorado River streamflow that have been implicated in the uncertainties of projections (Vano et al. 2014). First, our experiments are conducted at the nominally high global 50-km spatial resolution so as to capture high-elevation sources of runoff production with greater fidelity than other general circulation model (GCM)-based studies. Second, the spatial resolution achieved is fully dynamical in which the atmosphere is fully coupled to the land surface. This differs from prior methods that used either regional climate models or statistical downscaling of spatially coarse-grid GCMs as input to offline hydrology models (Gao et al. 2011; Christensen et al. 2004; Christensen and Lettenmaier 2007; Cayan et al. 2010; Vano et al. 2012). These offline uncoupled hydrologic models, though powerful tools for downscaling where topographic complexity is critical, may nonetheless introduce artifacts into hydrologic impact assessments from unresolved biases in downscaling and also from a distorted surface energy balance in LSMs that affects evapotranspiration (Pierce et al. 2013). Third, in light of considerable differences among LSM responses to identical meteorological changes, we address the fidelity of the global models’ various land model components. Following Milly et al. (2005) who studied future streamflow changes in global models of the Third Assessment of the Intergovernmental Panel on Climate Change (IPCC), we appraise in section 3 the climatological UCRB streamflow in our models. Additionally, a Budyko approach (Budyko 1974) is applied to appraise the realism of simulated water and energy balances that control linkages between runoff production and aridity in our models. The latter diagnosis serves to illuminate biases that could skew the response of streamflow to changes in aridity. Last, the atmosphere in our experiments is subjected to specified observed sea surface temperature and sea ice conditions [Atmosphere Model Intercomparison Project (AMIP) protocol]. Although not fully coupled to the ocean as in Climate Model Intercomparison Project (CMIP) simulations, our experiments incorporate realistic historical ocean states that are expected to reduce atmospheric biases relative to CMIP.
We present in section 4 the simulated UCRB hydroclimate responses to historical observed climate change driving. Using a non-fingerprint-based detection and attribution method (Wuebbles et al. 2017) that utilizes our ensemble modeling approach, the observed temperature, precipitation, and streamflow changes over the UCRB since the early twentieth century are interpreted in the context of climate change and variability. The streamflow change is further diagnosed in section 5 within a framework of its sensitivities to separate effects of surface air temperature and precipitation changes. A summary and concluding remarks appear in section 6.
2. Datasets and methods
Colorado River naturalized flow is an estimate of streamflow Q that would exist in a natural state unimpaired by diversions and withdrawals from human activities. We use annual water-year (1 October–30 September) estimates of such naturalized (“virgin”) flow at Lees Ferry based on Upper Colorado River Commission data spanning from 1896 to the present (see http://www.ucrcommission.com/RepDoc/UCRCAnnualReports/69_UCRC_Annual_Report.pdf). For the same period, water-year estimates of surface air temperature T and precipitation P are based on monthly 5-km gridded analyses (Vose et al. 2014), coarsened to a 50-km grid for comparison with our climate models. Although it is for a slightly shorter period (from 1906 to the present), we also use Bureau of Reclamation naturalized flow data (see https://www.usbr.gov/lc/region/g4000/NaturalFlow/Final-MethodsCmptgNatFlow.pdf), which includes information for eight subbasin catchments in the upper basin above Lees Ferry.
Changes in naturalized streamflow, temperature, and precipitation are calculated as differences between 30-yr averages for 1981–2010 relative to 1896–1925, the values of which are −13%, +1.1°C, and +1.5%, respectively.2 Uncertainty in estimating precipitation change is particularly large owing to the strong topographic dependency of precipitation and the large temporal and spatial inhomogeneities in observing stations (Milly and Dunne 2002), and this is evident in a comparison among gridded datasets (Henn et al. 2018). Sampling uncertainty is also of concern and will be addressed in section 5.
b. Climate models and experiments
Global atmospheric GCM experiments are diagnosed to determine UCRB sensitivity to historical changes in climate drivers. Three models are used including the National Center for Atmospheric Research Community Atmospheric Model (CAM5; Neale et al. 2012), the European Centre for Medium-Range Weather Forecasts/Hamburg model (ECHAM5; Roeckner et al. 2003), and Japan’s Meteorological Research Institute model (MRI3.2; Mizuta et al. 2017). Each represents comparable horizontal scales globally (~50 km over the UCRB), whereas they differ substantially in representation of land surface processes (see Table 1). Briefly, the CAM5 land model [Community Land Model, version 4 (CLM4)] is most sophisticated among the three, explicitly representing biogeophysical processes including surface radiation interactions with vegetation, stomatal physiology, and photosynthesis (Lawrence et al. 2011). The CLM4 includes 10 hydrologically active soil layers having a uniform depth of 3.8 m for active hydrology. In ECHAM5, a Simple Biosphere Model (SiB; Sellers et al. 1986) is used for calculating energy, mass, and momentum exchange between vegetated land surface and atmosphere. Its soil hydrology is further simplified to comprise a single layer of fixed 2-m depth, yet it allows for subgrid variability in soil infiltration capacity (Schulz et al. 2001). In MRI3.2, a SiB formulation is also used although involving three layers of active soil hydrology, with the lowermost layer set to 10 m (Hirai et al. 2007). We use the total of the models’ surface runoff and subsurface drainage term, that is, all terms that represent water leaving a model grid cell into the stream network, which on climatological time scales is nearly equivalent to streamflow. The CAM5 and ECHAM5 runs were performed by the authors, and the MIR3.2 runs were performed by the Japanese Meteorological Research Institute.
In the factual experiments, each model is forced by observed monthly SST and sea ice concentration variations using the boundary data of Hurrell et al. (2008) in CAM5 and ECHAM5, and Centennial In Situ Observation-Based Estimates of the Variability of SST and Marine Meteorological Variables, version 2 (COBE-2; Hirahara et al. 2014), in MRI3.2. Greenhouse gases (GHGs), aerosol, and ozone variability in CAM5 and MRI3.2 use the corresponding CMIP5 coupled model forcing protocols through 2005 (see Neale et al. 2012; Mizuta et al. 2017). For ECHAM5, GHGs vary according to the observed concentrations (Meinshausen et al. 2011), while tropospheric and stratospheric ozone vary based on Cionni et al. (2011). Aerosol concentrations do not vary interannually in ECHAM5. The CAM5 and ECHAM5 extensions after 2005 assume RCP6.0 forcing, while MRI assumes RCP8.5 forcing. The common period of these historical runs is January 1979–December 2010, and an ensemble of simulations is generated whereby each member of a particular model experiences identical time evolving boundary forcings but is initialized from different 1 January 1979 atmospheric initial states (see Table 1).
The counterfactual experiments are parallel runs in which the boundary conditions and atmospheric composition are specified as if global warming had not occurred since roughly preindustrial times. The procedures follow those used in extreme event attribution studies (e.g., Pall et al. 2011; Stott et al. 2012; Christidis et al. 2013, Massey et al. 2015; Sun et al. 2018). The models are forced with monthly varying boundary conditions that retain the interannual and decadal variability that occur in the factual experiment but in which estimates of long-term ocean warming have been removed. The GHG and ozone concentrations are set to late-nineteenth or early-twentieth-century values. Here, two different techniques are used in constructing the “pre–global warming” ocean boundary states—details of the model configurations and counterfactual development for CAM5 and ECHAM5 are provided in Sun et al. (2018), and those for MRI3.2 are in Mizuta et al. (2017). Briefly, for CAM5 and ECHAM5, counterfactual SSTs of 1979–2010 are generated by removing an observed 1880–2011 annual linear SST trend from the full variability. Only zonally averaged values of the trends are removed. For MRI3.2, a similar approach is used in that the SST warming pattern is also derived from observations except that the baseline for detrending is the early twentieth century, and trends at each grid point are removed rather than zonal averages. These two SST warming patterns are shown in Fig. 2 of Sun et al. (2018). Their common features are warming across all zonally averaged latitude bands having magnitudes of ~0.5°C. The regional structure of SST trends is characterized by greater warming in the Indian Ocean and tropical western Pacific relative to the tropical eastern Pacific (see also Solomon and Newman 2012). Note that removing a zonally averaged SST trend does not alter the zonal SST gradient present in the actual SSTs, whereas removing the full-field SST trends can. Given evidence from previous modeling studies for a considerable sensitivity of western United States precipitation to changes in the zonal SST gradient over the tropical Pacific (e.g., Hoerling et al. 2010), the manner in which the counterfactual SSTs are constructed is expected to affect the character of Colorado River basin precipitation. It is unclear what the best approach is for determining the pattern of global SST warming, and thus the best procedure for counterfactual construction. Herein, our use of two different approaches for generating boundary forcings absent global warming will accommodate some of the uncertainty. While noting similarity in large-scale features of meteorological changes across our experiments suggesting reproducibility in patterns (see the online supplemental material) across our three models, we have not treated the sensitivity of any single model to different counterfactual constructs and thus do not completely address the robustness of our results to various counterfactual approximations (see Christidis et al. 2013).
c. Diagnosing climate change impacts on UCRB hydroclimate
For 1981–2010, we calculate the factual minus counterfactual differences of UCRB temperature, precipitation, and runoff to quantify the effects of changes in climate drivers that have taken place over roughly the last century. These are compared with observed hydroclimate differences for 1981–2010 minus 1896–1925.
The ensemble modeling approach permits construction of a large sample of such centennial changes. Because the ensemble members are effectively independent samples of atmospheric variability, any combination of factual and counterfactual could occur with equal likelihood. For CAM5, a total of 100 (10 × 10) unique centennial differences are generated, whereas a total of 900 (30 × 30) and 10 000 (100 × 100) differences are calculated for ECHAM5 and MRI3.2, respectively, where the numbers in parentheses denote the number of factual and counterfactual experiments. The ensemble mean value of differences estimates the forced (signal) component of change while the spread among ensemble members estimates the internal variability (noise) component. We apply a simple detection and attribution approach to explaining the character of observed temperature, precipitation, and streamflow changes over the last century. The observed changes are compared with the statistics of simulated changes to determine the extent to which the former is consistent with variability that includes or excludes climate change drivers. This approach is referred to as a nonfingerprint-based approach to detection and attribution (Wuebbles et al. 2017). Following Wuebbles et al. (2017), an observed centennial-scale change is judged to have been detected when its magnitude has less than a 10% probability of occurring due to internal variability alone. This represents a detection of the signal with respect to atmospheric variability, with the oceanic (SST and sea ice) variability the same in each run, and as such underestimates the true natural variability. For brevity, the models are subsequently referred to as CAM, ECHAM, and MRI.
3. Simulated hydroclimates of the UCRB
Figure 2 compares UCRB observed water-year (1 October–30 September) temperature (Fig. 2, first column of top row) and precipitation (Fig. 2, bottom column of first row) with that simulated in the factual experiments of each GCM for 1981–2010. The principal common feature is a topographic organization of temperature and precipitation consisting of cold and wet conditions over the high-elevation eastern sections. By comparison, the western and southern lower elevations are relatively warm and dry. Despite these realistic spatial patterns in climatological conditions—owing largely to the models’ topographic fidelity—various biases are also apparent. Most noteworthy is each model’s 15%–20% wet bias above the observed basin-averaged water-year precipitation (Table 2). Each model also has a cold bias, ranging from only a few tenths of a degree Celsius for the CAM and ECHAM to nearly 2.5°C for MRI.
Climatological streamflow in each model (computed from the sum of surface runoff and subsurface drainage contributions as described in section 2) is also strongly controlled by orography in a manner broadly consistent with observations (Fig. 3). CAM and ECHAM generate high flow (nearly 60% of observed total basin flow) in the Yampa–White, Colorado Headwaters, and Gunnison basins consistent with observations (first column of top row in Fig. 3, dark blue shades). By contrast, and somewhat surprisingly given its otherwise realistic precipitation simulation, MRI (Fig. 3, top-right panel) fails to simulate appreciable flow in several high-elevation eastern catchments. The models agree overall with observations in having low streamflow in the western (Upper and Lower Green River and Dirty Devil) catchments. For the basin as a whole, excess streamflow beyond that observed at Lees Ferry is generated in all three models, qualitatively consistent with their positive precipitation biases. Specifically, as summarized in Table 2, the basin-average annual model streamflow ranges from 67 mm [totaling 15.3 million acre feet (hereinafter maf; 1 acre ft = 1233.5 m3) for the basin] in CAM to 80 mm (18.4 maf) in ECHAM. These compare to an observed Lees Ferry climatological flow of 64 mm (14.5 maf).
An additional metric of a basin’s hydroclimate is its streamflow efficiency defined as the ratio of streamflow to precipitation (Q/P). The observed streamflow efficiency over the UCRB as a whole is 16%, which compares to 14%, 18%, and 15% in CAM, ECHAM, and MRI, respectively (Fig. 3, bottom). Indicated hereby is the models’ realism in capturing a highly inefficient streamflow production over the basin whereby ~85% of total basin precipitation is lost to evapotranspiration. However, only CAM realistically captures the spatial heterogeneity in streamflow efficiency, being highest in the cold/wet eastern subbasins and lowest in the warm/dry western subbasins. The other models deviate appreciably from this pattern. For instance, while ECHAM has realistic streamflow efficiency in the eastern high-elevation subbasins, the efficiency is equally (and erroneously) high over the warm/dry western lower-elevation subbasins. The MRI streamflow fractions in the Colorado Headwaters and Gunnison basins (~10%) are notably unrealistic when compared with observations (~25%), despite having realistic precipitation magnitudes in those areas. This suggests an important mediating role played by land surface physics that also determine key features of streamflow. In addition to these climatological mean statistics, we have also compared times series of annual Lees Ferry naturalized flows to those simulated in the models. We confirm that the magnitude of the simulated interannual streamflow variability in individual model realizations is comparable to the observed variability during 1981–2010 (figure not shown).
We further assess the fidelity of GCM land surface physics using a Budyko framework that quantifies linkages between water and energy balances at catchment scales (e.g., Milly 1994; Sankarasubramanian and Vogel 2002; Milly et al. 2018). Figure 4 shows the broad features of a basin’s hydroclimate schematically by plotting the relationship between aridity [the ratio of potential evapotranspiration to precipitation (PET/P)] against runoff production [the ratio of actual evapotranspiration to precipitation (AET/P), equal to 1—Q/P where Q is total runoff]. Catchments that are characterized by large topographic diversity, such as the UCRB, can be populated by the full spectrum of hydroclimate regimes identified in the schematic. Small runoff production is typically associated with high values of PET/P that is characteristic of more arid regions toward the upper-right-hand side of the graph. In contrast, large runoff production is associated with low values of PET/P that is characteristic of wet regions (both regimes are illustrated in Fig. 4 by open circles). The Budyko framework applied at model gridcell scale across the UCRB thus yields a powerful process-based indication of hydrologic realism. However, there is not sufficient streamflow data at such fine scales to make observational comparisons with models feasible over the UCRB. Instead we compare GCM results with comparable high spatial resolution outputs from a hydrologic model [the Variable Infiltration Capacity (VIC) model; Liang et al. 1994] that has been driven with observed meteorological inputs for the UCRB, akin to using a reanalysis of streamflow (see the appendix). For the historical period of 1915–2015, the VIC-simulated annual UCRB runoff correlates with the naturalized Lees Ferry flow at 0.94, having a mean flow of 14.3 maf as compared with the Lees Ferry mean of 14.7 maf (see Fig. A2).
Figure 5 presents the Budyko analysis applied to each GCM. Starting at the left, CAM (red dots) shows a similar range to the observationally driven VIC model distribution affirming the vast contrast in hydroclimates across the UCRB. The results demonstrate the fidelity of physical linkages between aridity intensity and streamflow in CAM. ECHAM shows a somewhat weaker relationship between aridity and flow, resulting likely from tributary catchment-scale biases seen in Fig. 3. The MRI results are notably nonphysical, with several “humid” grid cells (red dots in the upper-left quadrant) generating almost no flow despite the model’s realistic temperature and precipitation. Overall, CAM simulates the intrabasin heterogeneity quite well, elements of which can also be discerned in the ECHAM simulations (Fig. 5, middle). By contrast, the actual evapotranspiration ratio in MRI is least realistic, differing substantially from the VIC results and also from those of the other two models.
Our working assumption is that models having a more realistic physical hydrology, as revealed through this Budyko diagnosis, may exhibit more realistic streamflow sensitivities to meteorological change (Milly and Dunne 2011). Guided by the vetting of the models as per diagnostics in Figs. 2–5, we will emphasize CAM hydrologic responses to climate change in subsequent sections. In contrast, vetting of our GCMs as concerns the likely realism of meteorological responses to global warming is less straightforward and less well informed by assessing climatological biases of UCRB basin temperature and precipitation. We note that contrary to the very different qualities of land surface models used in the three GCMs and their distinctly different biases, their atmospheric components have similar attributes including high spatial resolutions and a comparable realism of their climatological temperature and precipitation over the UCRB. None of these latter features necessarily dictate or discriminate how each model will respond to climate change forcing. In so far as this response may involve sensitivity to SST changes, we have verified that the GCMs exhibit comparably realistic atmospheric responses to the naturally occurring El Niño phenomenon (see also Zhang et al. 2016, 2018). Subtleties of global atmospheric biases that may bear upon the realism of a model’s regional Colorado River basin responses to climate change cannot be readily determined, however. We therefore will carry each of the three model’s temperature and precipitation responses to global warming forward in the subsequent evaluation because these will aid in informing the uncertainty in overall UCRB hydroclimate responses.
4. Simulated hydroclimate response of the UCRB
Figure 6 compares the statistics of UCRB averaged temperature (top), precipitation (middle), and streamflow (bottom) for 30-yr means (1981–2010) of factual experiments (red curves) to those of counterfactual experiments (blue curves). Differences in model statistics, for example as revealed by shifts in histograms, indicate the magnitude of responses to long-term climate change. As summarized also in Table 3, the effect of such climate change forcing is to cause temperature rise, precipitation decrease, and streamflow decline for the UCRB in each of our three models.
An important finding is that the detectability of such signals—indicated by the magnitude of forced change (differences between factual and counterfactual histograms) relative to internal variability (spread of the histograms)—differs appreciably among the three variables. For temperature, a complete separation of simulated 30-yr averages with and without climate change forcing is apparent for each model (Fig. 6, top). The warming signals range from 1.2° to 1.4°C and these are tenfold larger than the spread in statistics of 30-yr averages (Table 3). The models’ close agreement with the observed 1.1°C warming is thus consistent with such high signal-to-noise ratios. Our interpretation of these results is that a climate change–induced warming signal over the last century is highly detectable in the UCRB with respect to atmospheric variability, and that the observational analysis has indeed detected such warming.
The signal-to-noise ratio for precipitation change is much smaller, on the order of 1 relative to 10 for temperature. The signals range from −1.3% to −4.1% among the three GCMs, magnitudes that are comparable to the internal variability of 30-yr averages (Table 3). Consistent with such modest signal-to noise ratio, the precipitation histograms overlap considerably, although they are statistically different from each other (>99% with a Kolmogorov–Smirnov test) for the large ensemble of ECHAM and MRI runs (Fig. 6, middle). Conversely, the histograms are statistically indistinguishable, however, for the smaller ensemble of CAM runs. Concerning interpretation of the observed precipitation change, the +1.6% increase resides well within the range of model internal variability thus suggesting low detectability. The detection of precipitation change is further complicated by various systematic observational shortcomings in estimating the true observed precipitation change over the last century, primarily due to a limited and irregularly distributed station network (e.g., J. Barsugli et al. 2019, unpublished manuscript). Even in the case that this could be known accurately, each of the models indicates that a forced signal could be readily masked by internal variability. We thus judge a climate change signal in UCRB precipitation since the late nineteenth century to be undetectable at this time.
The streamflow signal (Fig. 6, bottom), although produced through a full energy and water balance at the surface, can be thought of as convolving each model’s land surface sensitivity principally with its temperature and precipitation changes. All models agree that meteorological changes in the Colorado River act to reduce flow. The magnitude of reduction is appreciable, and is close to the observed decline of −13% when averaged across the three models with a range in signals from −6.9% to −17.8%. The signal-to-noise ratio lies between 1 and 2 for streamflow, as compared to about unity for precipitation (see Table 3), and the streamflow histograms are statistically different from each other for each GCM. This indicates greater detectability for streamflow compared to precipitation due to the contribution of the highly detectable warming signal to runoff decline. We thus judge a 13% decline in UCRB flow over the last century to be marginally detectable (at a 90% level), and that the observational measurements of streamflow decline are unlikely to be due to natural variability alone. The notion of detectability during the instrumental record of a true change in UCRB streamflow must be tempered, however, by the fact that our model-based estimates of noise are only with respect to atmospheric variability, and thus conservative. Given the considerable decadal variability in flow volumes measured during the last century, and also from indications for strong decadal variability in the aforementioned nineteenth-century flow reconstructions, confidence of a detected change is judged to be low.
To provide finer granularity of simulated hydroclimate changes, we present in Fig. 7 the subbasin-scale responses for each GCM. Shown are the ensemble mean differences of factual minus counterfactual experiments (i.e., their forced signals). The pattern of temperature change is homogenous over the basin (Fig. 7, top) consistent with a large-scale warming pattern that encompasses the entire continental United States (see the online supplemental material). The pattern of precipitation changes (Fig. 7, middle) is more structured, consisting of a north–south gradient having largest declines over the southern basin. This dryness is part of a regionally focused precipitation reduction covering much of the U.S. Southwest that intrudes into the UCBR in each of the three models (see the online supplemental material). The pattern of streamflow changes (Fig. 7, bottom) shows declines in all subbasins (Fig. 7, bottom) and in all models, although having a more uniform pattern in ECHAM and MRI as compared with the more spatially variable pattern of CAM.
Several factors can account for the differences in streamflow responses among the three models. One arises from diverse land surface sensitivities of each model to temperature and precipitation. A second originates from the differences in meteorological forcing over the UCRB generated by each model. A few qualitative insights on these effects can be gleaned from Fig. 7, whereas a quantitative evaluation is provided in section 5. Notably, the larger flow declines in ECHAM and MRI simulations are plausibly due to their larger precipitation declines, relative to those in CAM. The magnitude of this effect will depend on the precipitation elasticity of streamflow, which could differ among the models. Applying an elasticity of ~2 based on observational data for this region (e.g., Sankarasubramanian and Vogel 2002, 2003), a −4.1% precipitation change signal in ECHAM versus −1.3% in CAM would imply streamflow change signals of −8.2% and −2.6%, respectively.
Streamflow sensitivity to temperature is also an important factor. Anticipating the more detailed analysis in section 5, one can already infer the role of warming by examining the CAM responses in Fig. 7. It is clear that its pattern of streamflow change does not resemble its pattern of temperature change, but instead more closely resembles its pattern of precipitation change (see Fig. 7, left column). If temperature were the primary driver of flow change, as some empirical studies propose, one would have expected the greatest streamflow decline in CAM to occur over the Upper Green basin where the greatest warming occurs. Instead, that catchment shows little streamflow change. The small signal of precipitation increases (~1%) in the Upper Green appears to be sufficient to offset effects of 1.5°C warming as concerns flow change. Implied here is a sensitivity to warming that is appreciably less than the values from −10% to −15% C−1 that are suggested from empirical studies. Alternatively, it is also possible that streamflow, at least in the Upper Green of CAM, is highly sensitive to precipitation—appreciably greater than an elasticity of 2. To clarify these issues, we therefore next diagnose runoff sensitivity to the separate effects of temperature and precipitation change.
5. UCRB streamflow sensitivity to temperature and precipitation change
To quantify the streamflow sensitivity to long-term meteorological change, we plot in Fig. 8 the scatter relationship between flow and precipitation percentage changes based on the factual minus counterfactual differences. The results indicate streamflow changes to be strongly constrained by precipitation changes in each GCM, despite various peculiarities of each model’s hydroclimates discussed in prior sections. Furthermore, their covariations are linear, having correlations exceeding 0.9. As such, the precipitation elasticity of streamflow can be reliably estimated by calculating the slope of the linear fit to the scatter (red lines in Fig. 8) the values of which are 2.4, 2.0, and 2.7 for CAM, ECHAM, and MRI, respectively (see Table 4). Such precipitation elasticity on centennial time scales is close to that derived from statistics of interannual variability of observed streamflow and precipitation (e.g., Sankarasubramanian and Vogel 2003) and is also close to that derived from idealized LSM experiments subjected to specified precipitation anomalies (e.g., Vano et al. 2012, 2014). The key point here is that precipitation change is critical for determining the long-term response of the Colorado River to climate change given this twofold amplification via elasticity.
The scatter relationships in Fig. 8 can be further diagnosed to also quantify the sensitivity of streamflow to temperature change. Note that the linear fit for each of the three scatterplots has a negative y intercept, indicating that flow is depleted even when precipitation change is zero. That depletion is reflective of long-term warming effects that occur in each model. When the y intercept is scaled for each model’s temperature change (see Table 3), the resulting temperature sensitivities of streamflow are −2.5%, −3.0%, and −6.5% °C−1 for CAM, ECHAM, and MRI, respectively. They are considerably lower than empirical estimates of UCRB streamflow sensitivity based on historical data analysis (e.g., McCabe et al. 2017) and also on the low end of LSM flow sensitivity to specified warming (e.g., Vano et al. 2012, 2014). For the roughly 1.1°C observed warming that has occurred over the last century, these high-resolution AMIP-derived sensitivities imply a temperature-induced streamflow decline of ~2.8%–7.2%. Note that these estimates of temperature sensitivity are for hydroclimate changes during the last century and may not be representative of future temperature–streamflow relationships.
Despite the high correlation of runoff and precipitation changes, the scatter around the linear fit is large when compared to the magnitude of the climate change signal itself. This scatter is a consequence of the different sequences of weather that are simulated in each ensemble member. Recognizing that observed estimates of runoff and precipitation change over the instrumental record constitute but a single point in that scatter diagram, there is considerable sampling uncertainty in knowing the true relationship between hydrologic and meteorological changes, even were each perfectly measured.
Figure 9 shows the spatial distribution of streamflow sensitivity based on analysis of the scatter relationships generated for each subbasin. Regarding temperature sensitivity (Fig. 9, top), reduced flow due to warming occurs across the entire UCRB in all models, although appreciable differences exist among the models in their spatial structures. For example, CAM shows a maximum sensitivity (~−4% °C−1) in the Colorado headwaters (eastern subbasins) where its maximum runoff production occurs. By contrast, ECHAM exhibits minimum temperature sensitivity (less than −2% °C−1) in those same areas. The MRI temperature sensitivity is also greatest in these high-elevation eastern subbasins, but with large magnitudes approaching −10% °C−1. Indeed, such large temperature sensitivities are mostly responsible for the overall high basin-averaged temperature sensitivity in MRI (see Table 4). Yet, we note that the MRI model was shown to produce erroneously little flow in these high mountain energy-limited regimes, contrary to observations (see Fig. 5). Owing to such severe climatological biases in its streamflow, the very high temperature sensitivity in those areas (and also its high temperature sensitivity for UCRB runoff overall) is judged to be unreliable.
For subbasin precipitation sensitivity (Fig. 9, bottom), CAM results indicate the high mountains of the Colorado River headwaters to be most sensitive to precipitation changes, with elasticities approaching 3. By contrast, neither ECHAM nor MRI models exhibit comparably strong spatial heterogeneity in their precipitation elasticities. Overall in CAM, these eastern catchments of high elasticity (>2.5) are locations of maximum climatological runoff production and are also the regions most sensitive to temperature change. Yet, owing to such locally high elasticity, even small precipitation changes would drive large flow responses that could easily overwhelm reductions in streamflow even under strong warming.
Multiplying these streamflow sensitivities with their corresponding temperature and precipitation changes (see Fig. 7), we present in Fig. 10 the individual contributions of temperature and precipitation to total flow change. The resulting linear additive reconstruction of total streamflow change (Fig. 10, top) closely reproduces the actual GCM changes (cf. Fig. 7 and Fig. 10) thereby affirming the merit of diagnosing each meteorological effect separately. The results demonstrate that diminished Upper Colorado River streamflow is primarily due to the effects of climate change–induced precipitation reductions. This is even true in CAM, which experienced the weakest precipitation reduction signal among the three GCMs, affirming the importance of precipitation elasticity. A second key result is that each of the GCM simulations indicates streamflow to be reduced by both a reduction in precipitation and a rise in temperature. Their combined effect thus leads to a substantially greater reduction in Colorado River flow volumes than each effect alone.
Last, we address the structural uncertainty in our assessment of the century-long decline in Colorado River flow and attempt to explain the appreciable range among the model signals (from −6.9% to −17.8%). Recall that the results of Fig. 10 intertwine two critical factors driving the Colorado River response to climate change over the last century, the first involving each model’s meteorological response to long-term climate forcing and the second involving each model’s land surface sensitivity to the associated meteorological changes. In aggregate, it is the manner in which these factors differ among the models that accounts for the structural uncertainties in streamflow responses. Concerning the first factor, we adopted a view that the meteorological changes in the three models have equal plausibility for reasons detailed in section 3. Concerning the second factor, section 3 did reveal distinct biases in water/energy balances that could bear upon land surface sensitivity. The Budyko-based diagnosis indicated CAM to be superior based on several climatological metrics, informing our view that its more realistic relationship between aridity and runoff production would yield a more realistic streamflow response to changes in aridity.
Given these considerations, we repeat the calculations of Fig. 10 but convolve the CAM sensitivity alone with each of the three GCM’s precipitation and temperature changes. The results of Fig. 11 can thus be viewed as comprising three estimates of streamflow change arising when the same land model (CAM) is subjected to three different scenarios of temperature and precipitation change. The principal features are not qualitatively different from those of Fig. 10 with basin-averaged streamflow declines of −6.4%, −13.9%, and −10.9% for CAM, ECHAM, and MRI meteorological changes, respectively.
6. Summary and concluding remarks
Colorado River flow has declined by roughly 20% during the instrumental record, as based on estimates of naturalized Lees Ferry flow. Using a high-resolution multimodel approach, we demonstrate that climate change forcing has acted to reduce streamflow from the Upper Colorado River basin during this period by about 10% (with uncertainty range of 6%–14% reductions). Our results thus indicate that about one-half of the observed flow decline during the instrumental period has likely resulted from long-term climate change.
A combination of factors stemming from the UCRB’s hydrologic sensitivity to changes in meteorological conditions account for how, and by how much, climate change has induced flow declines. Each of three different global models used herein indicate that climate change forcing during the last century has increased surface temperature (~+1.2°C) and decreased precipitation (~−3%) over the UCRB. The effect of precipitation change is shown to have been the most consequential factor owing to an amplified impact on flow resulting from precipitation elasticity of ~2. Warming has also driven runoff declines, although as a secondary factor over the past century owing to a modest −2.5% °C−1 temperature sensitivity based on our best model-derived estimate. We thus find that only about one-third of the overall climate change signal in Colorado River flow decline has resulted from warming, whereas about two-thirds of the climate change signal has resulted from precipitation decline.
It is useful to reconcile our findings with several recent studies that examined causes for low Colorado River flow since 2000, the so-called Millennium drought period. Udall and Overpeck (2017) used empirical data to estimate that about one-third of the anomalously low flow during 2000–14 (corresponding to a 19% deficit relative to twentieth-century climatology) was due to high temperatures. Based on their analyzed surface temperature anomaly of 0.9°C (see their Table 2), a −7% °C−1 runoff sensitivity can be empirically inferred. In a different approach, Xiao et al. (2018) used VIC simulations forced by an updated version of Hamlet and Lettenmaier (2005) meteorological analyses for 1916–2014 from which they determined that slightly more than one-half of the 2000–14 low flow was due to rising temperatures. Since their LSM simulations yielded only a 14% decline in Colorado River flow (relative to their model’s baseline historical simulation), a −7% °C−1 runoff sensitivity also emerges from their results. In both studies, the authors found a substantial fraction (from one-half to two-thirds) of low Colorado River flow during 2000–14 was due to precipitation deficits, the cause for which was subsequently found to be mainly internal atmospheric variability associated with natural fluctuations of tropical Pacific sea surface temperatures (Lehner et al. 2018). To synthesize, this collection of recent work indicates a substantial portion of the Millennium low flow regime in the UCRB has resulted from natural variability, but with a human-induced warming acting to also significantly deplete the river. Our own results, although not strictly addressing this drought event, are in qualitative agreement with a narrative that warming has exacerbated the low flow conditions since 2000.
The question remains open, however, with regard to the magnitude of that warming effect. The VIC model sensitivity to warming found in Xiao et al. (2018) is similar to that reported for a different set of VIC experiments examined by Vano et al. (2014). That study included intercomparison of eight different LSMs subjected to idealized UCRB warming the results of which indicated sensitivities at Lees Ferry ranging from about −3% to −8% °C−1. One of the stated purposes of our work was to overcome various limitations in such offline LSM experiments, and to this end we used coupled land–atmosphere GCM experiments to provide a more realistic treatment of surface energy and water balance adjustments to climate change forcing. Our best estimate of a −2.5% °C−1 temperature sensitivity is on the low end of these LSM results and is likewise appreciably lower than Xiao et al. (2018) and Udall and Overpeck (2017) based on their studies of the Millennium drought. Regarding the latter empirical study, sampling uncertainty can be an important limitation in determining the temperature sensitivity from observations alone. For instance, the correlation between annual temperature and UCRB flow is only about −0.5 (e.g., McCabe et al. 2017). Consistent with such a modest statistical linkage, our own examination of the Millennium drought that includes data thru 2018 (see time series of Fig. 1) indicates a more substantial fraction of the Lees Ferry flow deficit to be attributable to precipitation deficits (~80%) while temperature may have contributed only about 20% to the flow deficit. The inclusion of four additional years to the Millennium drought diagnosis reduces the estimated temperature sensitivity to −3% °C−1. This is not to argue that the latter is more representative of the true sensitivity, but merely to illustrate the considerable uncertainty in empirical approaches, some reasons for which are further pursued in J. Barsugli et al. (2019, unpublished manuscript). A more comprehensive study on the causes for the UCRB Millennium drought is clearly required, and because all trends calculated over the last century end in the Millennium drought event, the results would also bear on interpretation of the long-term flow decline.
We began our paper with the proposition that attribution would be tantamount to skillful prediction if Colorado River flow decline over the instrumental period had resulted principally from surface warming. To test that premise, we undertook a comprehensive examination of the past with a goal to better informing expectations of the future. It is our conclusion that temperature information alone may offer rather limited predictive value for future Colorado River flow, a consequence of our evidence for a weak temperature sensitivity derived thru analysis of a vetted multimodel large ensemble. For instance, only a −5% streamflow diminution would occur under a scenario of 2°C additional warming by the middle of the twenty-first century. Our findings show that such an impact could be readily masked by the inherent noise of internal variability. The message of our analysis is very clear—precipitation variations of the magnitudes that can arise naturally on multidecadal time scales would easily overwhelm such a warming effect for the foreseeable future given the large precipitation elasticity (>2).
A consequence of our finding is that the great matter confronting efforts to anticipate future Colorado River flow is how precipitation will change rather than how temperature will change. To be sure, the three models used herein were unanimous in indicating a drying effect of century-long climate change to date. But whether such an effect will prevail and grow as warming accelerates is unclear. We note that the median of CMIP5 model projections for RCP4.5 and RCP8.5 are for increases in annual precipitation over Colorado by the middle of the twenty-first century (Lukas et al. 2014; Ayers et al. 2016). If those projections come to pass, then Colorado River flow would plausibly increase, despite the further strong warming that is also anticipated. Several aspects of the precipitation change problem require further study, however. For instance, one key aspect toward predicting UCRB precipitation on such long time-scales involves the pattern of sea surface temperature warming in the tropical Pacific. The change pattern to date has deviated significantly from that simulated in CMIP models (Hoerling et al. 2010), and a science challenge is to understand why and what implications this has for reliability of existing projections over the UCRB. It should also be noted that our estimate of the basin’s runoff sensitivity to temperature rise was for warming over the last century—it is conceivable that the more substantial warming expected over coming decades may evoke a different sensitivity of basin runoff.
The study has addressed several sources of uncertainty in understanding how Colorado River streamflow responds to climate change [see Vano et al. (2014) for a list of many such uncertainties]. Our global model experiments were conducted at a nominally high spatial resolution (~50 km) so as to improve the representation of how topography organizes and controls the basin’s hydroclimate. Note that Vano et al. systematically investigated the effects of increased resolution on hydrologic sensitivity and found that higher spatial resolution results in weaker temperature sensitivity. Therefore, the relatively low temperature sensitivity of streamflow seen in our models is not likely due to course scale. A more robust appraisal of the realism of land surface physics was applied than in prior studies, helping to expose strengths and weaknesses of our models’ streamflow sensitivity to meteorological change. We employed fully coupled atmosphere–land modeling methods to avoid some of the biases that are introduced via techniques of statistical downscaling and the execution of offline hydrologic models. Also, our focus was on historical streamflow change rather than projections for the late twenty-first century (Christensen et al. 2004; Milly et al. 2005; Christensen and Lettenmaier 2007). Our experiments were therefore subjected to known oceanic and radiative forcing changes in the instrumental record during which the Colorado River was observed to decline. This distinction from studies that relied on model data derived from CMIP experiments is that the patterns of sea surface temperature change in those models appreciably differ from observations, leading to significant regional biases in North American climate and their sensitivity to climate change (Barsugli et al. 2006; Hoerling et al. 2010; Shin and Sardeshmukh 2011). However, using known changes in SSTs and GHGs (rather than projected values) allows for greater confidence in our assessment relative to future looking studies.
While having accommodated sources of uncertainty and addressed some of the ambiguities arising from model-based assessments, a full understanding of the century-long decline in Colorado River flow awaits treatment of various additional issues. One concerns the sensitivity of snowpack to climate change, and its role in determining the character of Colorado River sensitivity. Snow is especially relevant for the seasonality of Colorado River flow, which, however, was not a topic of this paper. Yet, seasonality in precipitation is believed relevant for water-year flow volumes given that cool season (October–April) precipitation explains 70% of the variability of UCRB water-year runoff since 1960 (Woodhouse et al. 2016). Observations indicate that snow water equivalent (SWE) on 1 April in the UCRB has declined during 1955–2014, and historical LSM runs indicate a longer-term SWE decline since the early twentieth century (Mote et al. 2018). It is unclear if these declines are symptomatic of a warming effect in the basin, and/or are due to precipitation declines—at least during the snow-accumulation season. It is important to emphasize that the true centennial-scale change in UCRB precipitation is not well known at this time (J. Barsugli et al. 2019, unpublished manuscript), and certainly not to the accuracy required for quantitative study of impacts on SWE and streamflow. Nor has the seasonality of its long-term changes been extensively explored, although this issue has been examined for decadal variability (Udall and Overpeck 2017). Continued research on the nature of historical observations including a reanalysis for long periods spanning the last century is thus important in its own right, but also toward better informing and more realistically constraining model investigations.
Land-use and land-cover changes over the United States as a whole contribute to local and regional changes in climate and water resources (Wuebbles et al. 2017). Perhaps most important is how vegetation cover and structure, forest growth, and land management practices such as crop versus wood harvest areas have changed in the UCRB over the last century. Neither of these factors are treated in this study where we used identical land surface/land-use conditions in our factual and counterfactual experiments. Bearup et al. (2014) and Livneh et al. (2015a) provide some evidence for hydrologic effects of forest transpiration loss on runoff that may have arisen in recent decades as a result of an infestation by bark beetles in the UCRB. However, these studies suggest that transpiration loss could serve to increase available water, which would not exacerbate the basinwide drying of interest. Additional evidence for land-use effects, although tied to remote land surface changes outside the UCRB, concerns dust on snow and the resulting consequences for Colorado River runoff efficiency, perhaps accounting for a 5% reduction in flows (Painter et al. 2007, 2010). A treatment of how land-use/land-cover changes over the entire UCRB have impacted Colorado River flows during the last century remains to be conducted.
Given these and other limitations, the current study cannot be viewed as comprehensive or as providing a definitive explanation for the decline in the Colorado River over the instrumental period. It does, however, offer a physically based explanation for the observed century-long decline in Colorado River streamflow within a new framework of the river’s sensitivity to meteorological change. We recognize that explaining this decline is urgent because Lakes Mead and Powell are now approaching levels (Barnett and Pierce 2008) that may imminently require delivery curtailments, as per the 2007 Interim Guideline Agreement for handling water shortages on the river (https://www.usbr.gov/lc/region/programs/strategies.html). To date, the impact of waning streamflow has been buffered by drawing down these reservoirs, although it is evident that the storage systems would become less of a bulwark against water shortage if the streamflow downturn persists. By providing this new approach toward assessing the sensitivity of the Colorado River to climate change, it is hoped that our results in combination with other studies will more robustly inform expectations about water resource abundance for coming decades and aid deliberations on interim guideline revisions, the current agreement of which expires in 2026.
We acknowledge Dr. Yukiko Imada for providing the data from the MRI3.2 climate simulations that were produced at the Meteorological Research Institute. We also acknowledge David Allured for generating the ECHAM5 climate simulations in the Physical Sciences Division. Comments on an earlier draft of this paper by Drs. Klaus Wolter and Robert Cifelli are also gratefully acknowledged. We also thank two anonymous reviewers and Dr. Dennis Lettenmaier for their comments on an earlier draft of the paper.
Land Surface Model Simulations
The VIC land surface model experiment is forced by gridded observed daily precipitation, maximum and minimum temperature, and wind speed based on Livneh et al. (2015b). Using soil parameters from Livneh et al. (2013, 2015b), the VIC model was built at a 0.5° resolution (~50 km) to conform with the AGCM scales, and we summarize its hydroclimate for the 1981–2010 period. The VIC simulation spans 1915–2015, and the historical experiments have been conducted at both ½° and 1/16° spatial resolution. The statistics of interannual variability and the correlation of the simulated historical annual flow with the natural Lees Ferry time series are very similar for these two resolution settings (Figs. A1 and A2).
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-19-0207.s1.
Early attempts at reconstructing Lees Ferry flow (LaRue et al. 1925; see also Meko et al. 2007) suggested that the decade of the 1880s may have had flow volumes as low as those in the 2000s decade, implying that the historical trend of Lees Ferry naturalized flow over a period of 1880–present may actually reveal less decline than during the instrumental record. Yet, LaRue’s reconstruction also paints a picture of large decadal variability in the nineteenth century, with indications that the decade of the 1860s may have been even wetter than before the compact.
Note that the streamflow decline of 13% computed in this manner is consistent with the linear trend that describes a 20% decline during 1896–2018 and a 24% decline during 1906–2018.