Observed snowpack loss

Observations of snow water equivalent (SWE) were acquired from the United States Natural Resource Conservation Service Snow Telemetry (SnoTel) network of automated snow pillow measurements across alpine sites. These measurements were taken from a quality-controlled data set discussed in ref. 12. We restrict our analysis to the post-1981 period since the number of observations before this time is too low to allow calculation of reliable regional averages. Monthly-mean (January-May) SnoTel observations are continuously available from 1982 to 2016 at 354 stations with elevations greater than 1,500 m. A threshold of 1,500 m provides good regional coverage; use of continuous observations avoids introducing temporal changes in spatial coverage. Figure 1 shows that 307 of the 354 stations (or about 87% of all stations) show a negative trend in annual maximum SWE (SWE max ). The maximum loss typically occurs in April. Here and subsequently, SWE max is computed over 1982–2010 to facilitate comparison with reanalysis-based estimates for the maximum period of overlap between the reanalyses and the SnoTel data (see below).

Figure 1: Topography and measurement network. The circles are the snow telemetry (SnoTel) network of stations utilized in this study. The red circles denote stations with negative trends in annual maximum snow water equivalent (SWE max ). The blue circles indicate stations with positive trends. Linear trends are computed from 1982 to 2010 to facilitate comparison with reanalysis-based estimates (for the period of maximum overlap between SnoTel and the four reanalysis products). All the stations plotted here are at elevations greater than 1,500 m. Full size image

We first compute the climatology and trend in SWE max . Results are area-averaged over the 354 selected SnoTel stations. Figure 2 shows that the climatology and trend are 4.6 cm and −0.33 cm per decade, respectively. The trend has a P value less than 0.15. This decrease in SWE max represents about a 9.5% loss; this is expressed as the per cent difference between the decadal averages of SWE max in 1982–1992 and 2000–2010 (the first and last decades of the period of maximum overlap). This multi-decade decline in regional snowpack is consistent with results from refs 13, 14, 15, despite the fact that these studies used different time periods, data sets and metrics. (The sensitivity to the metric used, for example, 1 April SWE or SWE max , has been discussed in ref. 16.) Finally, we note that the accumulated SWE over the winter months is primarily controlled by snowfall. Significant losses in SWE are observed throughout the snowfall season (Supplementary Fig. 1) and are not restricted to April, the month typically associated with SWE max .

Figure 2: Climatology and trend in area-average annual maximum snow water equivalent. The green bars show snow telemetry (SnoTel) ranges (±1 s.d.) computed from random samples of stations from the SnoTel network (see main text). The pink bars show ranges (±1 s.d.) from four land surface reanalyses (REANAL4). The black ellipse encompasses 95% of the climatology and trend values across 50 simulations of a global climate model. The circles represent the average climatology and trend for SnoTel, REANAL4 and the CanESM2 ALL simulations. All values in this figure are area averages of annual maximum snow water equivalent (SWE max ) at elevations greater than 1,500 m and north of 30°N and south of 50°N. The period considered is the maximally overlapping period from 1982 to 2010. Full size image

We also consider a group of four gridded SWE data sets obtained from reanalysis products (hereafter referred to as REANAL4). These four data sets have been assessed in ref. 17 over the common period of availability from 1982 to 2010. Whereas the SnoTel measurements provide a purely observational perspective, the REANAL4 group has been produced using snow schemes of varying complexity embedded in land surface models forced by meteorological reanalyses. The following data sets were used: (1) the Global Land Data Assimilation System Version 2 (ref. 18); (2) the European Centre for Medium-Range Forecasts Interim Land Reanalysis19; (3) the Modern Era Retrospective Analysis for Research and Applications20; and (4) the Crocus snow scheme driven by ERA-Interim meteorology21. Differences between these data sets have been discussed at length in ref. 17. The SWE fields were re-gridded to a 1° longitude by 1° latitude grid; monthly mean SWE results were averaged across the four reanalyses.

Area averaging SWE max at elevations greater than 1,500 m (north of 30°N and south of 50°N) yields a reanalysis-mean climatology and trend of 3.9 cm and −0.56 cm per decade, respectively. The trend has a P value less than 0.01. This decrease in SWE max represents an ∼21.8% loss between 1982–1992 and 2000–2010. The correlation between the SWE max time series from REANAL4 and SnoTel is 0.73. The range of climatology and trend estimates from the four reanalysis products encompasses the corresponding SnoTel estimates (Fig. 2). While this is encouraging, there are significant scale challenges in relating local SWE measurements to coarsely gridded snow reanalyses, and direct comparisons between the two must be approached cautiously22. The green bars in Fig. 2 show the uncertainties associated with the SnoTel values. Here, uncertainties are computed as ±1 s.d. across 10,000 random samples of the SnoTel stations. In each case, a station is randomly selected and the average climatology and trend is computed across a random selection of no more than 10 of its neighbouring stations within 1° longitude and 1° latitude (the spatial resolution of the reanalyses). In short, the SnoTel and REANAL4 results are consistent within the uncertainties in each data source; this holds for both climatologies and trends.

Historical simulated snowpack loss

Our ensemble of model simulations was generated with the Canadian Earth System Model version 2 (CanESM2; see Methods and ref. 23). Simulations include historical forcings from 1950 to 2005 and RCP8.5 forcing extensions after 2005. The RCP8.5 pathway is a high-emission scenario leading to an 8.5 W m−2 increase in radiative forcing by 2,100. Over the early part of the twenty-first century, forcing differences between RCP8.5 and other commonly analysed RCPs are relatively small24. Two 50-member ensembles were run with CanESM2, one with natural and anthropogenic forcings (ALL) and one with only natural forcings (NAT). Anthropogenic forcings include changes in well-mixed greenhouse gases, anthropogenic aerosols, tropospheric and stratospheric ozone, and land use. Natural forcings consist of changes in solar irradiance and volcanic aerosol loadings.

For the ALL ensemble, area-averaging the simulated SWE max at elevations greater than 1,500 m (north of 30°N and south of 50°N) yields an ensemble-mean climatology and trend of 6.1 cm and −0.57 cm per decade. The trend has a P value less than 0.001. The simulated trend is in good agreement with the average of the four reanalyses; the simulated climatology is at the high end of the reanalysis estimates (Fig. 2). The latter suggests greater climatological snowfall in the model than in reality, although this is difficult to assess given the absence of a reliable long-term observational record of winter precipitation. The simulated pattern of climatological SWE max is in reasonable agreement with the comparable field in reanalyses (see Fig. 3), albeit with spatial detail given the coarser resolution of the global coupled model (which is nominally 2.8° longitude by 2.8° latitude).

Figure 3: Climatology of annual maximum snow water equivalent. (a) Average of four land surface reanalyses. (b) Average of 50 global climate model simulations from the ALL ensemble. The period considered is that of maximal overlap between the reanalysis products (1982–2010). The orange contours are topographic height in kilometers. The outer contour is 1,500 m and the inner contours are in increments of 500 m. Full size image

Figure 4a shows 5-year mean anomalies of area-averaged SWE max at elevations greater than 1,500 m. Results are from SnoTel (green), REANAL4 (pink) and the ALL simulations (black). We also show results obtained by dynamical downscaling of 35 members of the CanESM2 ALL ensemble and RCP8.5 continuation (dashed curve). Downscaling relied on the Canadian Regional Climate Model version 4 (CanRCM4), which has a nominal resolution of 50 km (see Methods and ref. 25). Downscaling was performed over the North American domain defined in the Coordinated Regional Climate Downscaling Experiment26. The similarity between the CanESM2 and CanRCM4 ensemble-mean time series (r=0.99) is evidence that the temporal variability of area-average SWE max from CanESM2 is relatively insensitive to an increase in spatial resolution.

Figure 4: Anomaly of annual maximum snow water equivalent and scaling factor β. (a) Anomaly in non-overlapping 5-year averages of annual maximum snow water equivalent (SWE max ). Solid black is ALL ensemble-mean and grey is 10–90% range (based on 50 ALL realizations). Solid blue line is NAT ensemble mean and dashed blue lines indicate the 10 and 90% values. The ALL and NAT curves are from CanESM2. The dashed black curve is from the 35-member set of CanRCM4 simulations with ALL forcing. Pink denotes the average of four reanalyses and green is the SnoTel observations. Anomalies are defined in Methods. (b) Scaling factor and 5–95% uncertainty range estimated from application of an optimal fingerprint method to SnoTel observations, reanalyses and CanESM2 simulation output (see Methods). Scaling factors greater than 0 and consistent with 1 indicate that a model-predicted SWE max signal has been detected in observations and attributed to the imposed forcing changes in the ALL or NAT simulations. Full size image

Because there is only one realization of internal variability in the real world, we do not expect the CanESM2 ensemble-mean SWE max (which is averaged over many different model realizations of internal variability) to closely follow the time evolution of SWE max in the SnoTel data and REANAL4. The key point here is that the multi-decadal changes in SWE max are reasonably similar in SnoTel, REANAL4, and the CanESM2 and CanRCM4 ensembles under ALL forcing. Under NAT forcing, however, the CanESM2 ensemble mean does not replicate the observed decline in SWE max since the 1980s.

Detection and attribution analysis

We use a standard optimal fingerprinting method to address the causes of these SWE max changes. This involves regressing the observations onto the simulated ALL and NAT responses. The resulting values of the scaling factor β provide information on the extent to which the model SWE max responses must be scaled to best reproduce the observed SWE max changes. The method also yields the 90% confidence intervals on β estimates (see Methods).

Figure 4b shows that with combined anthropogenic and natural forcings, CanESM2 reproduces the reanalysis and SnoTel SWE max changes, albeit with smaller magnitude—that is, the scaling factors are significantly larger than 0 and less than (but still consistent with) 1. With natural forcings alone, however, the model does not reproduce the reanalysis and SnoTel changes. A similar result was obtained in ref. 13 using snow course data, observationally based precipitation, two earlier-generation climate models, two statistical downscaling approaches, and a late-twentieth-century analysis period (1950–1999).

Since we do not have a CanESM2 ensemble with anthropogenic forcing only, we compute the anthropogenic scaling factor by linear regression of the reanalyses or observations onto the difference between the ALL and NAT responses (denoted by ANT). We assume additivity of the ANT and NAT SWE max responses. Figure 4b shows that the estimated ANT SWE max changes: (1) are consistent with the reanalyses; (2) are detectable in the SnoTel observations but with significantly smaller magnitude; and (3) yield very similar β values to the ALL case. Regarding (2), when average SWE max is computed over the larger region encompassing elevations greater than 1 km (in order to reduce the noise contribution inherent in point measurements) anthropogenic forcing more successfully reproduces the magnitude of the SnoTel changes (Supplementary Fig. 2). While the results from this analysis provide clear evidence of an anthropogenic influence on snowpack water storage, it is difficult to more reliably quantify the magnitude of this influence given the relatively short observational SWE max record, the effects of underlying variability and the uncertainties inherent in our indirect estimate of the ANT SWE max response.

Near-term projected snowpack loss

Figure 5a shows 5-year mean SWE max from the ALL ensemble, averaged at elevations greater than 1,500 m. The first 5-year period is centred on 1988 and the last 5-year period is centred on 2038. The ALL ensemble uses the RCP8.5 forcing extension after 2005. The ensemble-mean response (black curve) shows a small increase in the rate of SWE max decline from the first half of the analysis period shown in Fig. 5a to the second half of the analysis period. More specifically, the externally forced trend in SWE max is about −0.50 cm per decade from 1988 to 2013 and about −0.62 cm per decade from 2013 to 2038. This increase in the rate of SWE max decline is significant at P<0.05, where significance is assessed using a standard difference of means test on the 50 simulated trends in the earlier period and the 50 simulated trends in the later period. As noted in ref. 15, distinguishing rates of change in externally forced response over small regions and short timescales, such as these, requires large ensembles of simulations (as shown here).

Figure 5: Annual maximum snow water equivalent and per cent change. (a) The black curve is for the CanESM2 ALL ensemble mean and shading is the 5–95% range. The ALL ensemble uses RCP8.5 forcing extensions after 2005. The red curve is the simulation with the greatest loss in annual maximum snow water equivalent (SWE max ) between 2013 and 2038. The blue curve is the simulation with the largest gain over the same 2013 to 2038 period. All values are 5-year averages plotted on the central year. (b) Per cent change in SWE max between the periods centred on 2013 and 2038 under different radiative forcing scenarios. Results in orange are from the CMIP5 multi-model ensemble (Supplementary Table 1). CMIP5 results are based on the analysis of one ensemble member (r1i1p1) from each CMIP5 model for which there exists a historical simulation with ALL forcing and a corresponding RCP2.6, RCP4.5 and/or RCP8.5 extension. The box-and-whiskers plots show ensemble-mean SWE max values, the 95% uncertainty ranges on the ensemble-mean values, and the minimum-to-maximum ranges. All values in this figure are for area averages of SWE max at elevations greater than 1,500 m and north of 30°N and south of 50°N. Full size image

Our focus is on quantifying the combined contribution of internal variability and external forcing to near-term projected SWE max . The largest near-term loss in accumulated snowpack over the western United States is ∼60%; this was calculated as the per cent change between averages over the 5-year periods centred on 2013 and 2038. The largest near-term gain is about 3%. The substantial difference between these two percentage changes—obtained with the same physical climate model, driven by the same external forcings—underscores the potentially large contribution of decadal variability to regional-scale changes in SWE max , particularly on shorter (ca. 30-year) timescales. The range in projected snowpack reductions reflects the interplay between temperature and precipitation trends on the end of season SWE15,27,28. Averaging over realizations reduces this variability, and provides a better estimate of the underlying response to external forcing, yielding about a 30% loss of regional snowpack in the ALL ensemble mean.

We note that CanESM2 generally reproduces the large-scale wintertime patterns of decadal variability observed in Pacific sea surface temperature, precipitation, sea level pressure and North American surface temperature11. This enhances our confidence in the credibility of the simulated decadal variations in snowpack loss. In comparison with other climate models participating in phase 5 of the Coupled Model Intercomparison Project phase 5 (CMIP5; Supplementary Table 1), the magnitude of unforced decadal variability in SWE max in CanESM2 is in the bottom half of the multi-model ensemble (Supplementary Fig. 3). This suggests that the CanESM2 estimate of the influence of decadal variability on regional snowpack is conservative. Alternately, smaller decadal variability also has implications for the detection results in Fig. 4, and may yield more liberal estimates of anthropogenic signal detectability.

We also compare the projections of regional SWE max changes in CMIP5 models and CanESM2 (Fig. 4b). Results are for changes in ensemble-mean simulated SWE max between the 5-year periods centred on 2013 and 2038; the minimum-to-maximum ranges are also indicated. The CMIP5 set consists of one ensemble member (r1i1p1) from each CMIP5 model with an ALL forcing historical simulation and corresponding RCP2.6, RCP4.5 and RCP8.5 extension (Supplementary Table 1). The range associated with the CanESM2 initial condition ensemble (for ALL forcing and the RCP8.5 extension) solely reflects the influence of decadal variability. The ranges associated with the CMIP5 multi-model ensemble are indicative of both decadal variability and model uncertainty.

Although the ‘spread’ in these ensembles arises from different reasons, the message from this comparison is that the potential for large snowpack loss in the CanESM2 initial-condition ensemble is also evident in the CMIP5 multi-model ensemble, and is manifest across all three emissions scenarios.

Large-scale patterns associated with extreme snowpack loss

Figure 6 shows the sea level pressure and surface air temperature patterns associated with a large negative contribution to regional snowpack from decadal variability in the period from 2030 to 2040 as obtained from the CanESM2 ensemble. These patterns are indicative of a stronger than normal Aleutian Low. The Aleutian Low is a semi-permanent low pressure centre located near the Aleutian Islands during the winter. It is one of the main centres of action in the atmospheric circulation of the Northern Hemisphere. The circulation anomalies associated with the stronger than normal Aleutian Low are responsible for transporting anomalously warm and moist air over the western United States. The impact of warm air on regional snowpack is obvious (especially where temperatures are above or about the freezing point). The impact of increased precipitation on snowpack is generally more complex29, but our model shows that the anomalous precipitation falls as rain rather than snow, therefore yielding a reduction in snowpack.