Dramatic declines in snowpack in the western US

Mountain snowpack stores a significant quantity of water in the western US, accumulating during the wet season and melting during the dry summers and supplying much of the water used for irrigated agriculture, and municipal and industrial uses. Updating our earlier work published in 2005, we find that with 14 additional years of data, over 90% of snow monitoring sites with long records across the western US now show declines, of which 33% are significant (vs. 5% expected by chance) and 2% are significant and positive (vs. 5% expected by chance). Declining trends are observed across all months, states, and climates, but are largest in spring, in the Pacific states, and in locations with mild winter climate. We corroborate and extend these observations using a gridded hydrology model, which also allows a robust estimate of total western snowpack and its decline. We find a large increase in the fraction of locations that posted decreasing trends, and averaged across the western US, the decline in average April 1 snow water equivalent since mid-century is roughly 15–30% or 25–50 km3, comparable in volume to the West’s largest man-made reservoir, Lake Mead. Mountain snowpack stores huge amounts of water in the western US, supplying much of the water used to grow crops. A team of researchers from Oregon State University and UCLA found that spring snowpack declined almost everywhere, especially in the coastal states and other locations with mild winter climate. (Skiers will be relieved that declines were smaller in winter.) Not surprisingly, the declines are mostly related to warming climate. Using a physically-based model of the hydrologic cycle, which takes daily weather as inputs and computes snow accumulation and melt, runoff, etc., the researchers computed the total snowpack in the western US. Total snowpack declined 15–30%, and the amount of that lost water is comparable in volume to the West’s largest man-made reservoir, Lake Mead. Many water managers are already planning for a future with less snow, but this research emphasizes that the future is here.


INTRODUCTION
California's recent multi-year drought (2011-16) and its extension into Oregon and Washington has shown that warming can create drought simply by preventing the accumulation of mountain snowpack. The year 2015, for instance, set the record low 1 April snow water equivalent (SWE) at over 80% of sites west of 117°l ongitude, 1 a result of high winter temperatures rather than low precipitation. [2][3][4] More than a decade ago, we showed that spring snowpack had declined at a large majority of locations in the mountainous western US, and corroborated the observations with hydrologic modeling that reached broadly similar conclusions. 5 We also noted that computing an area-averaged snowpack value from observations is challenging because the locations of long-term monitoring sites are usually chosen to favor a certain type of terrain and elevational range, with temperature-sensitive locations undersampled early in the record in some states. 6 Methodological choices (e.g., about record length) can therefore strongly influence results and must be carefully evaluated. In contrast, model-based estimates provide a basis for estimating long-term SWE changes across the entire Western U.S. domain.
Since our earlier work, several papers have further explored the relationships between mountain snowpack, variability and trends in precipitation and temperature, and geographically important factors. Stoelinga et al. (ref. 7 ) derived a snowpack index for the Cascades from streamflow measurements, from which they estimated that the spring snowpack declined 23% between 1930 and 2007. Pierce et al. (ref. 8 ) using a hydrologic model forced by observations and by two 1600-year climate model runs to estimate natural internal climate variability, attributed declines in snowpack (specifically SWE divided by accumulation-season precipitation) across the western US to anthropogenic warming. Luce et al. (ref. 9 ) postulated that decreases in westerly winds aloft may have contributed to orographically induced decreases in high-elevation precipitation and snowpack in the Northwest, which the sparse high-elevation observations might not detect. Our previous work found little to no decrease at high elevations, but large decreases at low elevations, 5,10 whereas their mechanism would produce the opposite. Nonetheless, the wind-induced changes join other factors besides direct anthropogenic warming in potentially influencing regional snowpack.
Changes in vegetation cover in or near snow courses, as well as station moves, could lead to spurious or non-climatic trends. Unfortunately, scant evidence is available about changes in vegetation over time at snow course locations. In a rare exception, Julander and Bricco (ref. 11 ) compare recent photographs of snow courses in Utah with photographs from 1936, noting that in some locations there was no change in vegetation while in others there was a complete transformation owing, e.g., to logging or fire. Such changes have been shown to affect significantly the accumulation and ablation rates, and could be as large as climate-driven trends. Julander and Clayton (ref. 12 ) in a study of 15 Utah snow courses (all but one at >2200 m elevation) noted that SWE tended to decrease where vegetation cover had increased. In the Supplementary Materials, we evaluate the effects of vegetation changes in the vicinity of snow courses in the Utah part of our domain using results of ref. 11 and show that while excluding stations that clearly have been affected by long-term vegetation change makes some difference to our inferences in specific locations, there is no clear relationship between changes in vegetation cover and changes in SWE, and the broad changes we find across larger portions of our domain are barely affected.
Precipitation and temperature can both play a role in determining year-to-year fluctuations in SWE and also in longterm trends. Mote et al. (ref. 5 ) showed that at the coldest sites, the correlations between winter temperature T and April 1 SWE are small and positive; trends there tend to be small. Moving toward warmer sites, the correlations between T and SWE generally decrease until they are large and negative; trends in SWE there tend also to be large and negative. They also noted that drier locations tended to be less sensitive to temperature. Luce et al. (ref. 9 ) mapped the relationships between SWE, T, and P explicitly by plotting April 1 SWE in a T-P space that spans western climates. They corroborated the sensitivity of SWE to T for warm locations found by ref. 5 and quantified the relationship of SWE to T and to P, providing a simple means to estimate how a given site's April 1 SWE would respond to a 3°C increase in temperature. Declines ranged from 9 to 100% (no snow) with a median 40%. Their analysis did not include the California Department of Water Resources data, so California locations are substantially underrepresented.
In this paper, we update our previous work using observed April 1 SWE data through 2016, and hydrologic modeling through 2014. As in Mote et al. (ref. 5 ) the temperature data used to drive the hydrologic model are adjusted to match the most reliable longterm stations (following the method outline in ref. 13 ) in an effort to ensure that long-term trends reflect best available knowledge of climatic trends. We also show results of an experiment using the hydrologic model in which warming trends were first removed from the model forcing data.

Spatial pattern of trends
Nearly all (92%) of snow courses and a large fraction (78%) of VIC grid cells now post negative trends over the updated period of record 1955-present (Fig. 1). In our previous work covering an earlier period of record, numerous sites in the cold, snowy, highelevation central and southern Sierra Nevada mountains in California had positive trends, in contrast to the VIC modeling which suggested only negative trends there. Now, California remains the state with the highest fraction of positive trends but extended drought in the most recent decade (post-2007) erased many of the formerly positive trends. Most states have only one or two locations with positive trends. Most of the largest negative trends are in eastern Oregon and northern Nevada, but trends < −70% also occur in California, Montana, Washington, Idaho, and Arizona. Only two sites have trends > + 70% and both have very low mean SWE (Big South in Colorado, 3.0 cm; and Siskiyou Summit, Oregon, mean SWE 3.8 cm). Conducting a Mann-Kendall significance test on the trends, we find that 232 of the 699 snow courses (33%) had significant negative trends at p < 0.05, and only 16 (2.2%) had significant and positive trends. For VIC, 17.5% of grid points had significant negative trends and 3.9% had significant and positive trends.
The observed and modeled patterns have similarities and differences, and there are also areas without observations to corroborate or contradict the modeled patterns. In Arizona, all of the observed and modeled trends are negative. There are no observations in the high terrain of the Mogollon rim area in northern Arizona or on the Arizona-New Mexico border just south of the Four Corners. In the San Gabriel and San Bernardino Mountains of southern California, VIC simulates declines at all grid points but there are no observations. In the observation-rich Sierra Nevada Mountains, there are pockets of upward trends, mostly at high elevations in the observations and mostly in the northern Sierras in VIC. In Oregon and Washington, declines predominate, but VIC suggests increases in the eastern Cascades of southern Oregon and the southern part of the Olympic Mountains (the two snow course locations are on the northern part of the Olympic Mountains). Trends in the Northern Rockies are almost uniformly negative in both model and observations, but they disagree strongly in eastern Idaho/western Wyoming where all observations show declines but VIC shows increases. In the Wasatch Mountains of Utah, VIC suggests strong increases, contrary to observations, but in Colorado, VIC uniformly shows decreases whereas a few observed snow courses show increases. Taken together, notwithstanding these local differences, in broad pattern both observations and model show large areas of decreases, with much smaller areas of increases.
In a simulation with VIC in which long-term trends were removed at each grid point from the model input temperature data ("no-warming"), long-term changes are dominated by changes in precipitation. Many areas with large negative trends in the regular VIC simulation (Fig. 1b) switch to smaller negative or even positive trends (Fig. 1c). Some of the largest changes occur in the Pacific states and northern Rockies, but most areas in Arizona also change from negative to positive trends. In other words, pointed out in the case of the Cascades, the mean elevation (and therefore sensitivity to temperature) of the network of snow courses decreased dramatically over the first few decades since early snow courses were mostly at higher elevation and the ones added later at lower elevations, which affects long-term comparisons and trends. This illustrates the importance of the criteria for selection of stations on the end results. Computing regional trends only from snow courses available before 1940 dramatically underestimates the temperature sensitivity and the magnitude of long-term trends, compared with snow courses available by 1950. California's snow course deployment was more evenly distributed in the vertical and from north to south, so the choice of starting year is somewhat less important there. We have not performed similar comparisons for the other regions. Although such distributional concerns do not apply to VIC, there are other reasons VIC outputs may not match the observations locally. Land use changes are not represented in VIC but there is some evidence that they may have affected long term trends (see discussion above and in the supplementary material about the work by Julander and Clayton (ref. 12 ).
Totaling the SWE across the VIC domain and converting to cubic kilometers provides a perspective on the magnitude of the changes (Fig. 3). The total decline of 21% in April 1 SWE is equivalent to 36 km 3 . For perspective, the capacity of the West's largest reservoir, Lake Mead, is 32 km 3 . Because the starting year of the VIC simulation and subsequent analysis is somewhat arbitrary, we compute linear trends of this time series in which the first year for the analysis (the "start year") is 1915, 1916, …1975 and the last year is always 2014. Although the slopes differ somewhat (Fig. 3b) the time series always has negative slope with declines that range from 15 to 30%.
We repeat that process for the other months. The seasonal pattern of trends is striking (Fig. 3b). January has the smallest trends for most start years. February and March are similar to April, though February has trends resembling January's for start years after about 1955. May and June have generally small trends for early start years but quite large downward trends (< − 30%) for start years after mid-century. For most months, the greatest percentage declines is for starting years in the late 1960s because the early 1970s were consistently snowy (Fig. 3a, Fig. 2).
These results for seasonal differences-generally becoming more negative from January through June -resemble those for the shorter satellite records of snow cover extent over the northern hemisphere, in which June had the largest declines of any month during the satellite era (1967-present). 15 Mapping each site and each model grid cell onto its seasonal mean temperature and precipitation, then averaging, the largest downward trends tend to be clustered in locations with mild wet climates (Fig. 4a) with mean temperatures above about −1°C and monthly mean winter precipitation >200 mm. Other climates where trends tend to be large and negative are found at locations with somewhat cooler drier climates (temperatures −5 to 0°C, precipitation <200 mm/mo). With VIC, similar patterns hold, with smallest trends roughly on a line from +3°C and 200 mm to −6°C and 300 mm and larger and more negative trends at both cooler  Fig. 2 Time series of regional mean 1 Apr SWE for the domains indicated in Fig. 1, for observations (black circles) and VIC (red crosses).
Smooth curves are added for VIC (red) and for the period of observations when at least a quarter (blue dashed) or half (blue solid) the locations were reporting. In computing regional averages for VIC, only grid cells with mean 1 April SWE > 50 mm are included drier climates and warmer wetter climates. VIC represents more low-elevation wet climates than the observations, and it is at these (temperatures >2°C and precipitation >400 mm/mo) where the largest negative trends are found, < −65%. Across the larger range of climates represented by the VIC domain for which there are no observations (outside the dashed box in Fig. 4b), the coldest locations tend to be quite dry and also to have predominantly large and positive trends, though for T < −15°C a number of the driest locations also have negative trends.

DISCUSSION
Since our 2005 paper, continued declines in western snowpack have resulted in an increase in the fraction of monitoring sites exhibiting negative trends. Independent calculations with the VIC hydrologic model broadly corroborate the decreases, and confirm that they are most prevalent in spring and are largest in the mildest locations. Differences between observed and VIC trends may result from (a) location and spatial distribution of observations (e.g., preference for siting snow courses on relatively level terrain and in key watersheds, lack of observations at low and very high elevations); (b) vegetation in the real world may change, whereas the VIC land cover is fixed through time in these simulations; (c) shortcomings in the VIC model forcing data (since station data are sparse at high elevations and generally in regions of complex terrain); and (d) errors in the VIC model's snow formulation (albeit these are most prominent during the ablation season, whereas we focus mostly on the snow accumulation season). Nonetheless, the spatial patterns, seasonal, and climatic dependence of the trends are broadly similar.
These declines are predominantly driven by warming trends, as confirmed by a VIC simulation with no warming trends, complementing earlier attribution work. 8 The declines in western snowpack represent a substantial loss of snow storage, as illustrated by our calculation in cubic kilometers showing that since 1915 western US snowpack has declined by 21% or 36 km, 3 greater than the volume of water stored in the West's largest reservoir, Lake Mead.
The magnitude of these changes relative to the built storage (reservoirs), and the certainty with which continued warming will lead to continued declines at a similar or increasing rate, 14 illustrates the immense challenge facing western water managers. Patterns of water use that became established (even entrenched) during the climate of the past cannot be changed without intense political effort owing to large cultural, economic, and infrastructure investments in the status quo ante. Solutions cannot consist solely of future infrastructure: new reservoirs cannot be built fast enough to offset the loss of snow storage, so solutions will have to lie primarily in the linked arenas of water policy (including reservoir operating policies) and demand management.

Snow observations
Our analysis uses data collected by the U.S. Department of Agriculture's (USDA's) Natural Resources Conservation Service (NRCS) (http://www.wcc. nrcs.usda.gov/snow/) and the California Department of Water Resources (CADWR) (http://cdec.water.ca.gov/snow/index.html) (see ref. 5 for details). These agencies operate hundreds of snowpack measurement sites -both manual and electronic -in the western US. Manual measurements are taken by snow surveyors along a "snow course", in which snow is measured by taking core and weighing the snow in the cores s (each reported "snow course" measurement typically is an average of multiple cores). Snowpack telemetry (SNOTEL) stations consist of a fluid-filled snow "pillow" and pressure transducer; SNOTEL stations automatically measure (and report) snowpack and weather data daily. Our work considers data from 1766 SNOTEL and snow course sites, screened first for longevity. A few snow courses have data going back to the early 20th century, including a few dozen in California that started before 1920, but most snow course observations began later in the century. The SNOTEL network was deployed beginning in 1979. In some cases, NRCS uses statistical relationships between SNOTEL and nearby snow courses to lengthen the records; such data values are marked in the data files as "estimated". In addition to NRCS and CADWR quality assurance procedures, we applied the following quality assurance checks: (1) check Apr 1/Mar 1 SWE ratios for validity (0.001 < ratio < 100), and (2) for the NRCS data (which included snow depth as well as SWE) the snow depth/ SWE ratio (a snow depth/SWE ratio < 1 was assumed to indicate invalid data). Only 18 data points across all years and all stations were found to be suspicious (and marked as such).
We use 1 Apr data because it is the most frequent observation date and early snow surveyors in turn selected 1 Apr because on average, the transition from the accumulation to the ablation season occurs around then, so it provides an estimate of the total water available for runoff. Data are nominally attributed to 1 Apr, but in reality, for some manual observations the closest measurement in a given year might have been collected some days before or occasionally after 1 Apr. CADWR data include an adjustment to 1 April, using nearby precipitation measurements and appropriate scaling relationships; we examined the trends in those data as well and found no significant differences compared with just using the date of observation. Stations selected for inclusion in this study had to have at least 50 years of data including at least 5 of the most recent years, and be at least 70% complete. In addition, we screened out some snow courses that we, as well as 12 Julander and Clayton (2015) identified as having significant land use changes (see Supplementary material). A total of 699 snow courses satisfied all our criteria.
For Fig. 2, the regional average shown in each panel was calculated for years when at least 20% of the eventually available long records were reporting. The smooth curve is shown as solid when at least 50% of the eventually available long records were reporting. These dates and the maximum number of snow courses were as follows: panel a, Cascades: 20% 1937, 50% 1952, maximum 128. b, California: 1930, 1939, 287. c, Rockies: 1937, 1951d, Dry interior, 1939d, Dry interior, , 1956, and 217.

Variable infiltration capacity (VIC) hydrologic modeling
We extracted SWE output from an updated version of the Hamlet and Lettenmaier data set 13 through 2014. VIC uses gridded daily precipitation and temperature (maxima and minima) data (primary variables), as well as secondary variables derived from the primary or other variables as described below. The extension to the data set made it current through 2014, and used the 1/16 degree latitude-longitude spatial resolution in contrast to the 1/8 degree spatial resolution in the original data set. The Hamlet and Lettenmaier approach uses a set of~1000 index stations from the Historical Climatic Network 16 across the conterminous U.S. which have lengthy records that have been corrected for the effects of instrument change and other effects that might otherwise result in spurious trends. These stations are augmented by a much larger set of NOAA Cooperative Observer stations with generally shorter records that help to define short term temporal and spatial variations. Long term variations (including trends) are adjusted to match HCN prior to gridding (details given in ref. 13 ) From the gridded precipitation and temperature data, we used algorithms described in ref. 17 to generate downward solar and longwave radiation, and surface relative humidity. Wind speeds are taken from the NCEP/NCAR reanalysis: 18 daily since 1949, monthly means before 1949. The primary and derived data are at 1/16 degree spatial resolution, to which the wind speed data are interpolated. We use these reconstructed, gridded data to force VIC, version 4.0.6, which produces SWE and other hydrologic variables. The VIC model is implemented using elevation bands to represent orographic variability in the VIC forcing variables. The number of elevation bands depends on the range of elevations within each grid cell; about 90 percent of the grid cells consist of a single elevation band; the remainder mostly have two bands, and a very small number have three.
The SWE values we use in our analysis are the average over the elevation bands.
Data availability