Testing the recent snow drought as an analog for climate warming sensitivity of Cascades snowpacks

We developed a 23 year climatology of SWE in the study watersheds for the period 1989–2011 using SnowModel, a spatially-distributed snow accumulation and energy balance ablation model (Liston and Elder 2006a, 2006b). We ran the model at 100 m horizontal grid spacing with a daily time step. To represent the distinct windward and leeward climates of the study watersheds, we defined separate modeling domains (figure 1) and defined unique temperature lapse rates for the MRB and UDRB. Lapse rates were computed from 30 arcsecond (800 m) PRISM monthly timeseries of minimum and maximum air temperature using standard least-squares linear regression, which were then interpolated between months and averaged to create daily lapse rates for each watershed. Precipitation (PPT) was partitioned into snowfall (SF) versus rainfall with a linear function of air temperature:

$$\begin{eqnarray}&&{\rm{S}}{\rm{F}}\\ =\,\left\{\begin{array}{ll}{\rm{P}}{\rm{P}}{\rm{T}},\, & {T}_{{\rm{a}}{\rm{i}}{\rm{r}}}\,\leqslant {\rm{T}}{\rm{S}},\\ \displaystyle \frac{1}{{\rm{T}}{\rm{R}}-{\rm{T}}{\rm{S}}}\times ({\rm{T}}{\rm{R}}-{T}_{{\rm{a}}{\rm{i}}{\rm{r}}})\times {\rm{P}}{\rm{P}}{\rm{T}},\, & {\rm{T}}{\rm{S}}\lt {T}_{{\rm{a}}{\rm{i}}{\rm{r}}}\lt {\rm{T}}{\rm{R}},\\ \,0, & {T}_{{\rm{a}}{\rm{i}}{\rm{r}}}\,\geqslant {\rm{T}}{\rm{R}},\end{array}\right.\end{eqnarray} \tag{ 1 }$$

where TS (−0.5 °C) is the temperature below which all precipitation is defined to be snowfall and TR (1.75 °C) is the temperature above which all precipitation is defined to be rainfall. The United States Army Corps of Engineers developed equation (1) from measurements of snowfall in the PNW (USACE 1956). We acknowledge the limitations of this empirical approach (Dai 2008, Marks et al 2013, Safeeq et al 2016) and provide a discussion of the uncertainty and how we handled it in the supplementary material.

Model inputs include (1) spatially distributed fields of precipitation, air temperature, wind speed, wind direction, and relative humidity; and (2) spatially distributed fields of topography and vegetation type. We used daily average air temperature, daily total precipitation, and daily average relative humidity gridded at 1/16° (5 × 7 km) (Livneh et al 2013), and the National Center for Environmental Prediction-National Center for Atmospheric Research (NCEP–NCAR) daily average wind speed and direction gridded at 2.5° spatial resolution (Kalnay et al 1996). Air temperature and precipitation were bias corrected with the PRISM 800 m monthly time series using the delta and ratio methods, respectively (Watanabe et al 2012). Wind and humidity data were bilinearly interpolated to the PRISM 800 m grid. The 800 m grids were used as the final model inputs, which were interpolated to the 100 m topography during model setup using SnowModel's built-in MicroMet utility. Topography and vegetation were resampled to 100 m resolution from the United States Department of Interior 30 m National Elevation Dataset digital elevation model and the 2011 National Land Cover Database land cover definitions (Gesch et al 2002, Fry et al 2011).

We used daily data from four automated SNOTEL stations ('test' sites) in the study region to evaluate SnowModel. The test sites span the windward and leeward sides of the study region and range in elevation from 1143 to 1733 m (figure 1). Some limited model development and calibration was required to achieve adequate model performance. This included adding a parameterization for retention and drainage of liquid water in the snowpack and a Monte Carlo parameter calibration routine (supplementary material).

Following Casola et al (2009), the effect of a change in surface air temperature on snow accumulation and melt can be assessed with a 'temperature sensitivity' parameter:

$$\begin{eqnarray}&&\delta x=\lambda (\delta T),\end{eqnarray} \tag{ 2 }$$

where $x$ is a snow-related metric, $\delta T$ is a change in surface air temperature, and

$$\begin{eqnarray}&&\lambda =\,\displaystyle \frac{{\rm{d}}x}{{\rm{d}}T}=\,\displaystyle \frac{\partial x}{\partial T}+\displaystyle \sum _{i}\displaystyle \frac{\partial x}{\partial {y}_{i}}\displaystyle \frac{{\rm{d}}{y}_{i}}{{\rm{d}}T}\,\end{eqnarray} \tag{ 3 }$$

is the temperature sensitivity of $x.$ The first term on the rhs of (3) is the 'direct' temperature sensitivity and represents changes in $x$ that result directly from changes in air temperature. The second term on the rhs of (3) represents indirect changes in $x$ due to a variable that is dependent on $T,$ for example an increase in local precipitation due to an increase in specific humidity with air temperature. In this study, we ignore all indirect effects and focus solely on the 'direct' temperature sensitivity, $\displaystyle \frac{\partial x}{\partial T},$ which we term $S.$ This choice simplifies the analysis for reasons of tractability, but also introduces uncertainty by ignoring potential indirect effects (e.g. Minder 2010). Following on this, if an estimate of $S$ is known, equation (2) can be used to predict $\delta x.$ For simplicity we adopt the notation:

$$\begin{eqnarray}&&D=S({\rm{d}}T),\end{eqnarray} \tag{ 4 }$$

where $D$ is a change in $x$ that results from a change in air temperature, ${\rm{d}}T,$ as predicted by an estimate of $S.$ In this study, we use a number of metrics for $x,$ described below.

We modeled $S$ using the 'delta' method (Hay et al 2000), where the historical daily air temperature data were uniformly increased by +2 °C (T2) and +4 °C (T4) and precipitation increased by +10% (P10). The temperature increases represent the expected mid- to late 21st century average changes for the Pacific Northwest (PNW) region from 20 global climate models and two greenhouse gas emissions scenarios (B1 and A1B) (Mote and Salathé 2010). For precipitation, the magnitude and direction of projected 21st century change for the PNW is highly uncertain (Safeeq et al 2016). Thus, a simple 10% precipitation increase was chosen to explore the effect of increased P on $S.$ The perturbed forcing data were used as input to the model to develop time series of SWE for each scenario, referred to as T2, T4, T2P10, and T4P10.

For each scenario we computed the mean per-pixel change in SWE (dSWE, m) and cumulative snowfall (dSF, m) on each calendar day, the mean per-pixel change in the date of peak SWE (dDPS, days), the snow disappearance date (dSDD, days), and duration of snow cover (dDSC, days), and the change in basin snow-covered area (dSCA, m²) on the date of peak basin volumetric SWE, relative to the reference period. These quantities were in turn used to calculate basin-scale area-weighted estimates of S for SWE, SF, and SCA on the basin DPS and 1st April, normalized by their reference period value and expressed as percent change (% °C⁻¹). For the timing metrics (DPS, SDD, SCD), we express estimates of S as absolute change e.g. dSDD/dT (days °C⁻¹).

To compare model findings with statewide SWE patterns during 2014 and 2015, we predicted statewide SWE anomalies with equation (4), using a model-derived empirical relationship between air temperature and S. We selected the metrics peak SWE, DPS, and SDD for this analysis. To achieve this, we binned grid cells by their mean DJFM (defined here as the snow accumulation season) air temperature, and computed the mean value of S for each 0.1 °C bin. We used these values as look-up tables that defined S as empirical functions of mean DJFM air temperature for each metric. We then computed the 2014 and 2015 PRISM DJFM air temperature anomalies to approximate ${\rm{d}}T$ and solved equation (4) for ${D}_{{\rm{S}}{\rm{W}}{\rm{E}}},$ ${{D}}_{{\rm{D}}{\rm{P}}{\rm{S}}},$ and ${D}_{{\rm{S}}{\rm{D}}{\rm{D}}}$ for PRISM grid cells with mean DJFM air temperature <4 °C and total DJFM precipitation >15 cm (values that bound the range of maximum DJFM air temperature and minimum 1 April cumulative precipitation computed from the SNOTEL data):

$$\begin{eqnarray}&&{D}^{2014}={S}^{{\rm{T}}2}\times {\rm{d}}{T}^{2014},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{D}^{2015}={S}^{{\rm{T}}4}\times {\rm{d}}{T}^{2015},\end{eqnarray} \tag{ 6 }$$

where $S$ is the look-up table sensitivity for each metric and ${\rm{d}}T$ is the PRISM DJFM air temperature anomaly. To evaluate the robustness of the extrapolation, predicted $D$ was compared to observed 2014/2015 peak SWE, DPS, and SDD anomalies $({D}_{{\rm{o}}{\rm{b}}{\rm{s}}})$ for each SNOTEL station in the state of Oregon with a period of record beginning 1981 or earlier.

Previous studies have used empirical methods e.g. linear regression between interannual variations in air temperature and SWE ('temporal analogs'; Mote 2006, Stoelinga et al 2010, Luce et al 2014), physically-based model simulations (Sproles et al 2013), or a combination of the two approaches (Mote et al 2005, Casola et al 2009) to estimate SWE temperature sensitivity in the PNW region. Luce et al (2014) presented a third method, termed 'spatial analogs' that used nonlinear polynomial regression between spatial variations in air temperature and SWE. These studies have led to some disagreement (Mote et al 2008, Casola et al 2009, Luce et al 2014), thus it is important to note that our approach is both similar (e.g. the physical model) and unique (e.g. the extrapolation procedure) from these previous methods. We discuss similarities and distinctions between ours and previous findings in section 5. Luce et al (2014) provide a detailed discussion of previous work on the topic.

Basin-wide modeled declines in peak SWE and the concurrent loss of SCA for the four scenarios are shown in figure 2. Peak basin-mean SWE decreased by 59% and 88%, and SCA on the basin DPS decreased by 40% and 69%, for the T2 and T4 scenarios, respectively. The 10% increase in precipitation had a relatively small effect on decreases in SWE and SCA, but the effect was elevation (and thus temperature-regime) dependent (figure 3). In our relatively warm basins, the basin-mean SWE losses were nominally offset by 7% and 2%, and SCA losses by 1% and 1%, for the T2P10 and T4P10 scenarios, respectively. Expressed as basin-mean sensitivities, ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}2}=29 \% ,$ ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}2{\rm{P}}10}=26 \% ,$ ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}4}=22 \% ,$ and ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}4{\rm{P}}10}=22 \% .$ Shifts in timing metrics (DPS, SDD, DSC) were also large and insensitive to increased precipitation (table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  The basin-mean sensitivity of each metric for each scenario, and the mean modeled sensitivity at basin test sites compared to observed sensitivities for 2014 and 2015. Snowfall is not measured at the test sites and snow-covered area is not a point-scale metric.

	Basin-mean				Test-sites
	${S}^{{\rm{T}}2}$	${S}^{{\rm{T}}2{\rm{P}}10}$	${S}^{{\rm{T}}4}$	${S}^{{\rm{T}}4{\rm{P}}10}$	${S}^{{\rm{T}}2}$	${S}^{2014}$	${S}^{{\rm{T}}4}$	${S}^{2015}$
Peak SWE (%)	−29	−26	−22	−22	−26	−33	−21	−22
1 April SWE (%)	−31	−27	−23	−22	−28	−35	−23	−25
Peak SF (%)	−28	−25	−21	−21	—	—	—	—
1 April SF (%)	−26	−24	−21	−20	—	—	—	—
Peak SCA (%)	−20	−19	−17	−17	—	—	—	—
1 April SCA (%)	−25	−23	−21	−20	—	—	—	—
DPS (days)	−13	−13	−13	−13	−11	2	−11	−19
DSC (days)	−34	−31	−31	−29	−28	−32	−24	−21
SDD (days)	−25	−23	−23	−22	−21	−32	−19	−22

These modeled changes in SWE were primarily controlled by shifts from snowfall to rainfall, both on a basin-scale (table 1), and across the entire elevation range (figures 3(a) and (b)). The largest relative decreases in SWE (figure 3(a)) occurred at low elevations where absolute decreases (figure 3(b)) were small. The largest absolute decreases occurred at elevations between 1700 and 2500 m where the majority of SWE accumulated during the reference period. Changes in cumulative snowmelt (figure 3(b)) appeared to be a consequence of the reduced SWE as opposed to a cause, as total snowmelt decreased at elevations up to 1700 m. Above 1700 m, snowmelt increased slightly, but increased melt was small relative to the reduction in snowfall. The strong dependence of these changes on modeled shifts from snowfall to rainfall is demonstrated by the basin-scale mean wet-day air temperature distribution (figure 3(d)). The nonlinearity of modeled $S$ suggested by comparing basin-scale ${S}^{{\rm{T}}2}$ to ${S}^{{\rm{T}}4}$ is at least in part explained by the step-function shift toward warmer temperatures we apply to this distribution, however the true effect will depend on changes to the statistical moments of this distribution in both space and time under future climate.

During 2014 and 2015, the mean DJFM air temperature anomalies at the four SNOTEL test sites were +1.3 °C and +3.9 °C, respectively. Basin-averaged PRISM air temperature anomalies were +1.0 °C and +3.5 °C, respectively. Cumulative SNOTEL precipitation for 2014 and 2015 on 1 April was 103% and 92% of normal, whereas basin-averaged PRISM precipitation was 97% and 90% of normal, respectively. These unique meteorological conditions roughly mirror what is expected for the region by mid- to late 21st century (Mote and Salathé 2010) and thus motivate our comparison between the modeled T2/T4 SWE, measured SWE during 2014/2015, and corresponding sensitivities (figure 4).

Zoom In Zoom Out Reset image size

During 2014, the peak SWE anomaly was −40%, DPS was 3 days later than normal (1 April), and SDD was 39 days earlier than normal (3 June). The T2 peak SWE anomaly was −53%, DPS was 24 days earlier (9 March), and SDD was 41 days earlier (1 June). During 2015, the peak SWE anomaly was −83%, DPS was 73 days earlier than normal (16 January), and SDD was 82 days earlier than normal (21 April). The T4 peak SWE anomaly was −85%, DPS was 39 days earlier (23 February), and SDD was 80 days earlier (23 April). Expressed as sensitivities, ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{2015}$ (−22%) was strikingly similar to our estimate for ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}4}$ (−21%), but ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{2014}$ was underestimated by ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}2}.$ The timing metrics ${S}_{{\rm{D}}{\rm{P}}{\rm{S}}},$ ${S}_{{\rm{S}}{\rm{C}}{\rm{D}}},$ and ${S}_{{\rm{S}}{\rm{D}}{\rm{D}}}$ showed varying levels of agreement between modeled and observed sensitivity but the modeled ${S}_{{\rm{D}}{\rm{P}}{\rm{S}}}$ was notably inconsistent with observations during both years (table 1).

Further examination of figure 4 shows that temporal SWE variability during 2014 and 2015 corresponded to the sub-seasonal phasing of precipitation and temperature anomalies (which the mean modeled sensitivities necessarily smooth out), shedding light on the deviation between observed and modeled ${S}_{{\rm{D}}{\rm{P}}{\rm{S}}},$ ${S}_{{\rm{S}}{\rm{C}}{\rm{D}}},$ and ${S}_{{\rm{S}}{\rm{D}}{\rm{D}}}.$ This behavior was consistent across sites during 2014, when early winter SWE was below average but peaked 3 days later than normal on 1 April, and during 2015 when the opposite phasing occurred, SWE peaked 73 days early on 16 January, and melted out 82 days early on 21 April, as opposed to more modest shifts predicted by ${S}_{{\rm{D}}{\rm{P}}{\rm{S}}}$ and ${S}_{{\rm{S}}{\rm{D}}{\rm{D}}}$ (table 1). The effect of storm phasing was particularly strong in 2014 when a very cold and wet February boosted record low January snowpacks and the resulting 1 April SWE was recorded as the 8th lowest in the 35 year period of record. Each site showed a pattern of recovery that resulted in a normal DPS. Contrasting behavior was observed during 2015 when early season precipitation coincided with anomalously warm temperature and the basin-wide temperature anomaly remained above 3.4 °C from January through March, with no recovery as in 2014. As a result, SWE peaked on 16 January and the lowest 1 April SWE in 35 years was recorded despite 90%-of-normal 1 April precipitation.

These results demonstrate that 2014/2015 peak SWE in our study watersheds closely resembled the T2/T4 peak SWE, but also that observed and modeled DPS and SDD were less consistent in both absolute terms and when expressed as S. Here, we show how the basin-derived S predict observed D across a broader geographic area during 2014 and 2015. We first note that figure 3 suggests $S$ strongly depends on elevation, a result noted by many (Knowles and Cayan 2004, Mote 2006, Tennant et al 2015). However, this dependence is largely due to the dependence of air temperature on elevation (Mote 2006). Thus, a more general approach is to define $S$ as a function of mean air temperature during the snow accumulation season. We use these relationships (figure 5), as opposed to the basin-mean sensitivities, to extrapolate our estimates of $S$ across the state (figure S2), and to test the degree to which observed conditions during 2014 and 2015 correspond to those predicted by the modeled sensitivities.

Zoom In Zoom Out Reset image size

Statewide, $\,{D}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{2015}$ predicted by ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}4}$ corresponded reasonably well to observed anomalies during 2015 (r2 = 0.52; figure 6(d)) when winter conditions were unusually warm (+3.3 °C at >1000 m) and precipitation relatively normal (−1% at >1000 m), with minimal spatial variability throughout the state (figure S3). During 2014, precipitation variability (−17% at >1000 m) likely contributed to the underpredicted peak SWE anomalies (r2 = 0.39; figure 6(a)). In general, there was much more variability in the observed response to conditions in 2014 and 2015 than the mean sensitivity functions predicted (figure 6). Only peak SWE showed similar spread across the observations as the extrapolated predictions. The timing metrics DPS and SDD were particularly scattered. These results are consistent with the plot scale analysis where the phasing of precipitation appeared to determine the temporal character of SWE, which likely differed from phasing across the state. Further, the 1σ spread in the sensitivity functions (figure 5) suggest substantial uncertainty across relatively narrow climatic bands, especially at warmer locations and for greater magnitudes of warming.

Zoom In Zoom Out Reset image size

When averaged across all sites, the 2014 statewide mean peak SWE anomaly was −36%, DPS was 9 days earlier than normal on 13 March, and SDD was 52 days earlier than normal on 19 April. Mean values predicted by ${S}^{{\rm{T}}2}$ were −30%, 17 days, and 27 days (for comparison see yellow squares, figure 6). The 2015 mean peak SWE anomaly was −67%, DPS was 56 days earlier than normal on 26 January, and SDD was 101 days earlier than normal on 1 March. Mean values predicted by ${S}^{{\rm{T}}4}$ were −68%, 40 days, and 64 days. The cool temperatures and late precipitation during spring 2014 resulted in a DPS later than predicted, but the SDD occurred 25 days earlier than predicted. The early precipitation and warmth during winter 2015 resulted in a DPS 26 days earlier than predicted, and the SDD 37 days earlier than predicted. Expressed as sensitivities for the Oregon SNOTEL network, ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}2}$ = 29% °C⁻¹, ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}\,}^{2014}$ = 30% °C⁻¹, ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}}^{{\rm{T}}4}$ = 20% °C⁻¹, and ${S}_{{\rm{S}}{\rm{W}}{\rm{E}}\,}^{2015}$ = 20% °C⁻¹.