Colorado River flow dwindles as warming-driven loss of reflective snow energizes evaporation

See allHide authors and affiliations

Science  13 Mar 2020:
Vol. 367, Issue 6483, pp. 1252-1255
DOI: 10.1126/science.aay9187

Evaporating futures

Drought and warming have been shrinking Colorado River flow for many years. Milly and Dunne used a hydrologic model and historical observations to show that this decrease is due mainly to increased evapotranspiration caused by a reduction of albedo from snow loss and the associated rise in the absorption of solar radiation (see the Perspective by Hobbins and Barsugli). This drying will be greater than the projected precipitation increases expected from climate warming, increasing the risk of severe water shortages in an already vulnerable region.

Science, this issue p. 1252; see also p. 1192


The sensitivity of river discharge to climate-system warming is highly uncertain, and the processes that govern river discharge are poorly understood, which impedes climate-change adaptation. A prominent exemplar is the Colorado River, where meteorological drought and warming are shrinking a water resource that supports more than 1 trillion dollars of economic activity per year. A Monte Carlo simulation with a radiation-aware hydrologic model resolves the longstanding, wide disparity in sensitivity estimates and reveals the controlling physical processes. We estimate that annual mean discharge has been decreasing by 9.3% per degree Celsius of warming because of increased evapotranspiration, mainly driven by snow loss and a consequent decrease in reflection of solar radiation. Projected precipitation increases likely will not suffice to fully counter the robust, thermodynamically induced drying. Thus, an increasing risk of severe water shortages is expected.

The Upper Colorado River Basin (UCRB) supplies water to ~40 million people and supports ~16 million jobs (1). Atmospheric warming and recent precipitation deficits have heightened concern about the future (26), but the response of river discharge to warming remains highly uncertain. An implicit assumption in the literature on UCRB hydroclimatic change is that two climatic mean variables—precipitation and temperature—determine runoff (hence, river discharge) response, following constant sensitivities α [percent discharge change per percent precipitation change (dimensionless)] and β [percent discharge change per degree Celsius of warming (% °C−1)]. Empirical regression analyses imply large values of β (−13 to −15% °C−1) (4, 68), which is inconsistent with estimates in the range −2 to −9% °C−1 obtained from perturbation of temperature inputs (the delta method) to hydrologic model simulations (2, 9, 10) and from theory (11). For α, regression and delta estimates are in much better agreement (10). The discrepancy in β, which is seen for rivers around the globe (11), translates into great uncertainty in the magnitude of future effects on human livelihood, economic activity, and ecosystem health. The situation is exacerbated by limited process understanding in the presence of hydroclimatic nonstationarity (12). The empiricism that is inherent in the regression approach, and even that which is inherent in the estimation of energy-driven evaporative demand in the hydrologic models (13), leaves the use of such methods for extrapolation of past observations to the future, under anthropogenic climate change, open to question. Accordingly, we gave special attention to surface net radiation—the ultimate driver of evapotranspiration—and to its modulation by snow-affected surface albedo (14) rather than relying on temperature measurements as a surrogate for energy availability. We found a strong influence of snow-affected albedo on radiation balance in the UCRB (Fig. 1) (15), which necessitated its consideration in a process-based estimation of β.

Fig. 1 Observed relations among monthly SWE, surface albedo, and surface net radiation in the UCRB.

(A) Dependence of surface albedo on SWE (logarithmic scale) for each of 12 elevation ranges. 1st, 2nd, and 3rd quartiles of binned data are shown. Curves are least-squares fits to the unbinned data and are used in the model. (B) Inferred dimensionless sensitivity C¯Rn¯dRndC of net radiation Rn to co-albedo (one minus albedo) C as a function of mean elevation of 960 subareas by month of the year. Blue curves are fitted to smoothed (30-point moving median; black) data from empirical regression estimates. Red curves are analogous fits for theoretical case where a change in absorbed solar radiation causes no radiative feedbacks. Fits to regressions are used in the model, except that fits to no-feedback data are used in a sensitivity experiment.

Herein, we address the following questions, in turn, by use of a monthly water-balance model grounded in a suite of observations: Does the model reproduce the historical regression-based β? What is the model’s delta-based β, and why does it differ from the regression-based value? Can the two values be reconciled? What physical processes control β? How sensitive is our β estimate to the assumptions in our analysis? How much did warming contribute to the historical hydrological drying in the UCRB? What future changes in UCRB discharge can be expected?

In addition to the snow-water equivalent (SWE), albedo, and radiation measurements used to develop the relations in Fig. 1, we used observations of precipitation and temperature (Fig. 2 and Fig. 3, A and B) as well as discharge (Fig. 3G) to constrain a hydrologic simulation model (15), which we used to elucidate the processes that control sensitivity and to reconcile divergent, previously published sensitivity estimates. The model has a monthly time step and divides the 290,000-km2 UCRB into 960 subareas to capture the strong heterogeneity induced by rugged (2700-m relief) topography (Fig. 2C). Rain-snow partitioning depends on temperature. Evaporative potential is set to the rate of non–water-stressed evapotranspiration under conditions of minimal advection (16). Fifteen model parameters were estimated by maximizing goodness of fit to observed discharge (15). We measured goodness of fit with respect to mean, linear trend, regression-based sensitivities α and β, and Nash-Sutcliffe coefficient of efficiency. [Including a correction that accounts for temporary subsurface storage of runoff before entering the river (11), which has previously been neglected, and using an October to September water year, we found observational regression-based α and β, ± one standard error of estimation, to be 1.98 ± 0.16 and −16.1 ± 2.9% °C−1, respectively. Neglecting the storage correction yields β = −13.1 ± 2.4% °C−1, consistent with earlier analyses.] The sensitivity of our results to the goodness-of-fit criteria is presented in the supplementary materials (15).

Fig. 2 Spatial distributions of annual climate variables and elevation over the UCRB, as resolved by the model discretization of space into 960 subareas.

(A) Total precipitation. (B) Mean temperature. (C) Elevation.

Fig. 3 Water-year time series of basin-mean, annual-mean values.

(A to G) Precipitation (millimeters per year) (A), temperature (degrees Celsius) (B), April 1 SWE (millimeters) (C), surface albedo (D), surface net radiation (watts per square meter) (E), evapotranspiration (millimeters per year) (F), and discharge per unit area (millimeters per year) (G). Blue curves represent estimates from observations, and gray bands represent ensemble range of model outputs. Least-squares linear fits also are shown.

Of 500,000 trial parameter sets, 171 satisfied the goodness-of-fit criteria (15), and these formed a model ensemble for subsequent analyses. As the temperature rose, the ensemble-mean SWE and—hence, following the relations in Fig. 1—albedo decreased, which led to a basin-mean increase of net radiation by 3.0% per century over the study period (Fig. 3, C to E). With an associated increase in evapotranspiration (Fig. 3F), the ensemble mean–modeled annual discharge (Fig. 3G) fell by 20.1% per century, compared with 19.6% per century observed; the square of the correlation coefficient (r2) between the observed and ensemble mean–modeled annual discharge is 0.82. Within the ensemble, the models’ regression-based, storage-corrected sensitivities α and β ranged from α = 1.89 ± 0.16 to 2.08 ± 0.18 (mean, 1.99) and β = −15.4 ± 2.9 to −16.9 ± 3.0% °C−1 (mean, −15.9% °C−1), which is consistent with observational estimates.

The ensemble was rerun with the temperature increased by 1°C, and differences from the base simulations were used to estimate sensitivities. The delta-based β ranged from −7.8 to −12.2% °C−1 (mean, −9.3% °C−1), which is consistent with higher-magnitude values from previous delta-based analyses (2, 10) and lower than regression-based estimates. Simulations with precipitation perturbed by +1% each month yielded a delta-based α of 2.21 to 2.83 (mean, 2.52). It is not surprising that some previous delta estimates of β were as high as ours nor that some were lower, because models had a variety of features, with varying physical realism, differentially sensitizing their evapotranspiration to temperature, and the mechanisms underlying β generally were neither identified nor constrained with measurements.

We found that the difference between the model’s regression- and delta-based sensitivities is explained by confounding variables: seasonal shifts of precipitation that historically accompanied annual anomalies of temperature. During warm water years, precipitation tended to shift from December–April to August–September. Because discharge sensitivity to precipitation is strong in winter (when extra precipitation tends to run off) (3) and weak in summer (when extra precipitation tends to evaporate), runoff is suppressed during warm years. To remove the confounding variables from our comparison of the delta and regression estimates, we modified the original experiments so that the monthly course of precipitation in every year was set to climatology times a factor that preserved the observed annual anomaly; for temperature, the annual anomalies were applied as additive constants to the climatology. For these experiments, the models’ regression-based α and β were α = 2.26 ± 0.03 to 3.00 ± 0.08 (mean, 2.59 ± 0.04) and β = −6.5 ± 0.7 to −11.6 ± 1.0% °C−1 (mean, −8.1 ± 0.7% °C−1), which is in reasonable agreement with the delta-based values; the difference between −8.1 ± 0.7 and −9.3 is only marginally significant, and allowances must be made for the simple formulation of the storage correction for the regression estimate (15).

The substantial dependence of inferred sensitivities on seasonal distributions of climate perturbations implies that the use of simple annual sensitivity parameters (α and/or β) can severely distort climate-change analyses. This is a shortcoming of both the regression and delta approaches. With regression, the derived sensitivities depend on basin-specific historical intra-annual and interannual variability, including the confounding precipitation-temperature covariance. In the usual delta approach, the perturbations have no seasonal variations, and the roles of precipitation and temperature are decoupled, so delta sensitivities are more readily interpreted. However, the best approach for hydroclimatic projections is to use the delta approach with projected monthly varying climate changes.

To understand the ensemble-mean magnitude −9.3% °C−1 of the delta-based β and its potential relevance for ongoing anthropogenic climate change, we consider the physical processes at play. Temperature enters the model in four ways: (i) Because SWE depends on the phase of precipitation and the rate of snow melt, the surface albedo and, hence, the evaporative potential are temperature-dependent. (ii) The maximum fraction of net radiation that is converted to latent heat flux depends on the temperature-dependent slope of the saturation vapor-pressure curve (11, 16). (iii) Evapotranspiration from soil ceases below a critical temperature, simulating winter dormancy of vegetation. (iv) Temperature affects the timing of snow melt and, thus, causes differences in sublimation and evapotranspiration in the model. By disabling these processes one at a time, we found that the contributions from the first three processes were −6.2, −2.1, and −0.3% °C−1, respectively, and other snow-storage effects and nonlinear interactions accounted for the remainder. Figure 4 summarizes the foregoing reconciliation of sensitivity estimates.

Fig. 4 Summary of estimates of β from previous studies and from this analysis.

Left to right: previous observational (Obs) and model analyses (2, 4, 610) and results from this analysis. Error bars represent ± one standard error of estimation from the regressions. The multicolored bar shows the contribution of each of the temperature-sensitive mechanisms to the magnitude of β. Excluded as unrealistic from the previous delta analyses are cases in which maximum daily temperature was perturbed whereas minimum was not (10). The label “P&T have anomalies only at annual scale” refers to the computations in which the monthly course of precipitation in every year was set to climatology times a factor that preserved the observed annual anomaly.

We repeated the analysis under the assumption that a change in albedo induces negligible radiation feedbacks (Fig. 1B, red). We found an ensemble mean β of −7.8% °C−1, indicating that our findings are somewhat sensitive to uncertainties in albedo-radiation feedback.

Unaccounted factors in our analysis include externally driven changes in radiation (e.g., from changing atmospheric composition), changes in boundary-layer entrainment (17), and stomatal responses to CO2 fertilization (18). The latter two factors tend to decrease the efficiency of the conversion of net radiation to potential evapotranspiration. We found the potential net effect of these factors on β to be negligible (15).

Our parameterization of potential evapotranspiration by use of the Priestley-Taylor formulation (16), which allowed for no atmospheric aridity feedback caused by actual (nonpotential) evapotranspiration, could be questioned. We therefore repeated our analysis with allowance for this feedback (15), finding a negligible difference in results. Another caveat to consider is that our adoption of the Priestley-Taylor formulation, even when we consider the aridity feedback, implicitly assumes that variabilities (in particular, long-term trends) of wind speed and humidity will not affect the value of the Priestley-Taylor α, even though they do, on certain time scales, play a documented role in variabilities of pan evaporation (19) and of the American Society of Civil Engineers Standardized Reference Evapotranspiration (20).

How much have temperature changes contributed to the period-of-record discharge trend (−19.6 and −20.1% per century observed and modeled, respectively) and the 2000 to 2017 discharge deficit (−15.9 and −17.6% of previous mean observed and modeled, respectively)? If we set temperature every year to its climatology, the model yields a discharge trend of −8.4% per century and a discharge deficit of −8.1%. We conclude that temperature sensitivity accounts for more than half of both drying phenomena, which is consistent with a previous analysis (21).

What about the future? To characterize future temperature and precipitation, we used the 8 out of 24 Coupled Model Intercomparison Project Phase 5 (CMIP5) climate models that simulated 1913 to 2017 discharge (area-weighted runoff) within a factor of two of the observed discharge. (The constraints used for climate-model selection were much less stringent than those used for selection of hydrologic model parameter sets because the hydrologic model was driven by historical climate time series, whereas the modeled climate time series were biased and independent of actual history.) From CMIP5’s historical, Representative Concentration Pathway 4.5 (RCP4.5), and RCP8.5 scenarios, we computed month-of-year temperature climatology increases from 1913–2017 to 2036–2065, added them to the observed historical record, and reran the ensemble. Across the set of eight climate models, ensemble-mean discharge decreased 14 to 26% (RCP4.5) and 19 to 31% (RCP8.5).

Could possible future increases in precipitation counteract the temperature-driven drying? When month-of-year temperature increases and precipitation ratios from the climate models were both applied, the ensemble-mean discharge decreased 5 to 24% (RCP4.5); under RCP8.5, changes ranged from an increase of 3% to a decrease of 40%. Thus, it appears unlikely that precipitation changes will be sufficient to fully counter the temperature-induced drying, though they might moderate it.

Many water-stressed regions around the world depend on runoff from seasonally snow-covered mountains, and more than one-sixth of the global population relies on seasonal snow and glaciers for water supply (22). It has been well established that snowpack serves as a reservoir that beneficially regulates the timing of water availability (23). Our findings imply that snow cover is also a protective shield that limits radiation absorption by, and consequently evaporative losses from, this natural reservoir; incidentally, this explains the observed phenomenon of precipitation as snowfall favoring runoff (24). The progressive diminution of this ecosystem service as a result of climate change will have a deleterious effect on water availability in snow-fed regions that are already stressed, including the UCRB.

Supplementary Materials

Materials and Methods

Supplementary Text

Fig. S1

Tables S1 to S5

References (2532)

Data S1 to S3

References and Notes

  1. Methods and supplementary text can be found in the supplementary materials.
Acknowledgments: This study was facilitated by the Geophysical Fluid Dynamics Laboratory of the National Oceanic and Atmospheric Administration and by several data providers cited in the supplementary materials. The authors gratefully acknowledge colleague reviews by R. Koster and T. Delworth. Funding: The authors are supported by the U.S. Geological Survey. Author contributions: P.C.D.M. was responsible for the conceptualization and overall direction of the work and wrote the original draft. P.C.D.M. and K.A.D. carried out computations. K.A.D. performed data curation and reviewed the original draft. Competing interests: The authors declare no competing interests. Data and materials availability: No original data collection was performed. The results of this study are reproducible and extensible by use of the cited data sources and other information in the supplementary materials.

Stay Connected to Science


Navigate This Article