This study examines the changing roles of temperature and precipitation on snowpack variability in the Northern Hemisphere for Second Generation Canadian Earth System Model (CanESM2) historical (1850–2005) and future (2006–2100) climate simulations. The strength of the linear relationship between monthly snow water equivalent (SWE) in January–April and precipitation P or temperature T predictors is found to be a sigmoidal function of the mean temperature over the snow season up to the indicated month. For P predictors, the strength of this relationship increases for colder snow seasons, whereas for T predictors it increases for warmer snow seasons. These behaviors are largely explained by the daily temperature percentiles below freezing during the snow accumulation period. It is found that there is a threshold temperature (−5±1°C, depending on month in the snow season and largely independent of emission scenario), representing a crossover point below which snow seasons are sufficiently cold that P is the primary driver of snowpack amount and above which T is the primary driver. This isotherm allows one to delineate the snow-climate regions and elevation zones in which snow-cover amounts are more vulnerable to a warming climate. As climate projections indicate that seasonal isotherms shift northward and toward higher elevations, regions where snowpack amount is mainly driven by precipitation recede, whereas temperature-sensitive snow-covered areas extend to higher latitudes and/or elevations, with resulting impacts on ecosystems and society.
Snow cover is a major component of the climate system that can extend over one-third of the global land area during boreal winter (Dozier et al. 1989). Because of the high albedo and low thermal conductivity of snow, as well as its influence on hydrology and biogeochemical cycles, receding snow cover due to climatic warming will have major impacts on climate feedbacks, ecosystems, and society (IPCC 2014).
Snow cover is a temporally and spatially integrated response to snowfall and snowmelt events (Brown and Mote 2009). Sufficiently cold air temperatures favor solid precipitation and snow accumulation, whereas downwelling shortwave and longwave radiation fluxes provide most of the energy enabling snowmelt (Sicart et al. 2006). Because warm air temperature plays a role in melt through turbulent heat exchange and is thus a good proxy indicator of snowmelt, snowpack variability can be largely inferred from near-surface air temperature and precipitation variability alone (e.g., Anderson 1973; Cayan 1996; Hamlet et al. 2005; Mote 2006; Brown and Mote 2009; Beniston 2012; Morán-Tejeda et al. 2013; Sospedra-Alfonso et al. 2015). Consequently, the large-scale response of snow cover to a warming climate is largely determined by the projected changes in air temperature T and precipitation P (IPCC 2014). At higher latitudes or elevations where cold season temperatures are well below freezing, snowpack amounts may not change or may increase as projected increases in snowfall can offset the shorter snow accumulation period. In lower latitudes or elevations where average winter temperatures are closer to the 0°C isotherm, additional warming results in large reductions in snowpack amounts from the combined influences of a shorter snow accumulation period, a decreasing fraction of solid precipitation, and more frequent winter melt events.
The sensitivity of snow cover to climate effects, particularly the relationship between snow and climatic variables, has been studied previously. Hantel et al. (2000) employed station observations to show that the relationship between snow-cover duration in Austria and the seasonal mean temperature over Europe adjusted with the elevation of the measurement stations can be described with a logistic function. This methodology was later implemented for Switzerland (Wielke et al. 2004) and improved and extended to Europe (Hantel and Wielke 2007). Morán-Tejeda et al. (2013) employed station measurements in the Swiss Alps to examine the changing influences of P and T on snow depth as a function of elevation. They highlighted the linear dependence between elevation and the potential for P and T to control snow depth interannual variability. Sospedra-Alfonso et al. (2015) carried out similar studies for the central Rocky Mountains and also found a linear dependence for the potential of P and T to control snow water equivalent (SWE) for the elevation range and time period investigated. More recently, Scalzitti et al. (2016) employed dynamical downscaling and climate projections to assess how a changing climate may impact the relative roles of P and T in controlling SWE variability in the western United States. Because of the high density and large elevation range in their model data, a logarithmic dependence was derived between elevation and the correlation between P predictors and SWE. Other studies on the roles of P and T in driving snowpack variability and trend include those by López-Moreno (2005) and Beniston et al. (2003a,b) for alpine regions in Europe; Cayan (1996), Hamlet et al. (2005), Knowles et al. (2006), and Nolin and Daly (2006) for alpine regions in North America; and Qin et al. (2006) and Ma et al. (2010) for the Tibetan Plateau. Although these studies help characterize temperature and precipitation shifts as indicators for changes in snow amount and/or duration, they have been limited to specific, mostly alpine regions, and, in some cases, such interpretation would implicitly assume that the climate and its natural variations are stationary.
To assess global-scale trends in a changing climate, Räisänen (2008) examined projections of SWE in 20 models from phase 3 of the Coupled Model Intercomparison Project (CMIP3) and showed that changes in precipitation, solid precipitation ratio, and snowmelt efficiency are explanatory of the simulated changes in SWE across the Northern Hemisphere. He emphasized the strong dependence of these projected SWE changes on present-day temperature and showed that the −20°C isotherm of November–March (NDJFM) mean temperature in the late twentieth century largely divides the snow cover into regions with increasing and decreasing midwinter SWE. Also, Räisänen (2008) showed that regions with NDJFM mean temperature below −25°C will very likely experience SWE increase in the late twenty-first century, whereas SWE is unlikely to increase in regions with NDJFM mean temperature above −10°C.
The above findings are consistent with those of Brown and Mote (2009), who summarized previous results describing trends in several snow quantities and found that the response of snow cover to a warming climate varies strongly with geographic location, elevation, time period, season, and snow-cover variable (e.g., snow-cover duration, snow depth, snow-covered area). Through a sensitivity analysis of a snowpack model and the study of observed and simulated changes of snow cover, they identified snow-climate regions and elevation zones that are expected to receive the earliest impacts among these diverse responses in a changing climate. Using the National Centers for Environmental Prediction (NCEP) reanalysis (Kalnay et al. 1996) in 1961–90 they showed that SWE has particularly high sensitivity to warming in moist regions where snow season temperatures range between −5° and +5°C. They also found a general consensus in the CMIP3 multimodel dataset for significant increase of SWE by 2070–99 in regions where November–May mean air temperature was below −20°C in 1970–99, whereas significant SWE decrease was found in regions having mean snow season air temperatures near the 0°C isotherm.
Although these results consistently find a strong dependence of SWE changes to present-day temperature, the functional relationship and quantification of such dependence is lacking. In this paper, we propose a framework to describe the geographic (e.g., latitude/longitude, elevation) and temporal dependence of SWE response to changes in P and T, using Second Generation Canadian Earth System Model (CanESM2) climate simulations and the National Aeronautics and Space Administration (NASA) Modern-Era Retrospective Analysis for Research and Applications (MERRA; Rienecker et al. 2011) as verifying observations. We derive a functional relationship between the influence of P or T on snowpack amount variability and the climatological seasonal temperature, defined here as the climatological mean near-surface air temperature during the snow accumulation period (a more precise definition is given in section 2e). Since the timing of snow onset and melt generally depend on geographic location and year, the term “seasonal” refers to time periods that are generally different in different regions and/or epochs. Consequently, this functional dependence is derived without the assumption of a stationary climate and is found to be largely independent of time and emission scenario.
The opposing dependence of snowpack response to P and T variability with mean seasonal temperature enables us to identify a climatological isotherm largely independent of time and emission scenario, which divides the snow cover into two distinct regions: one where snowpack amounts are mostly precipitation driven and another where snowpack amounts are mostly driven by temperature. In turn, this result enables us to delineate regions and elevation zones where snowpack accumulation is more vulnerable to a warming climate and show that, as climatological isotherms shift northward, there will likely be a significant reduction of precipitation-driven snow-cover extent by the end of the present century. Our approach differs from previous efforts in that the analytical expressions derived for the changing influences of P and T on snowpack variability apply to the entire Northern Hemisphere and are irrespective of emission scenario. These expressions can thus be adapted to more regional studies, including those previously cited examining dependences on elevation in mountainous regions, as well as specific climatic conditions.
The remainder of this paper is organized as follows. Section 2a describes CanESM2, MERRA, and the data employed. Section 2b defines the snow and climatic indices, whereas the role of climatic indices as predictors of the snow indices is described in section 2c. A multiple linear regression model and measures of the degree to which P and T drive snowpack variability are described in section 2d. The use of seasonal temperatures as proxies of geographic locations is described in section 2e. A basis for quantifying the relative effect of P and T on snowpack as a function of seasonal temperature and for identifying a threshold isotherm dividing snow cover into two distinct regions having different snowpack behaviors is given in section 2f. Results are presented and discussed in section 3, and a summary and conclusions are provided in section 4.
2. Data and methods
CanESM2, developed at the Canadian Centre for Climate Modelling and Analysis (CCCma), combines the physical components of the Fourth Generation Canadian Coupled Global Climate Model (CanCM4; Merryfield et al. 2013) with an interactive carbon cycle (Arora et al. 2011). CanESM2 provided CCCma’s long-term climate simulations for phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012) that informed the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC).
The atmospheric component of CanESM2 is the Fourth Generation Canadian Atmospheric General Circulation Model (CanAM4), described by von Salzen et al. (2013). CanAM4 is a T63 spectral model with physical fields represented on a ~2.8° Gaussian grid and 35 vertical levels with a lid height of ~1 hPa. The ocean component of CanESM2 combines the Canadian Model of Ocean Carbon (CMOC; Zahariev et al. 2008) with the Fourth Generation Canadian Ocean Model (CanOM4) described by Merryfield et al. (2013). The land component is version 2.7 of the Canadian Land Surface Scheme (CLASS2.7; Verseghy 2000) with land carbon described by version 1.0 of the Canadian Terrestrial Ecosystem Model (CTEM; Arora et al. 2009; Arora and Boer 2010). Sea ice is simulated using the cavitating fluid model of Flato and Hibler (1992). Radiative forcing includes the influence of greenhouse gases and aerosol emissions, explosive volcanoes, and solar variability as described by Arora et al. (2011).
We employ data from five CanESM2 ensemble members for daily SWE, near-surface air temperature, and precipitation. The data are analyzed in a 30-yr window moving in 1-yr increments over 1850–2100, encompassing the CMIP5 historical and RCP 2.6, 4.5, and 8.5 scenarios. Following Brown and Mote (2009), we define a snow year extending from 1 August to 31 July and exclude grid cells having SWE mm prior to 1 August to avoid snow persistence between years (i.e., we only consider seasonal snow cover). In addition, we exclude regions having monthly SWE interannual standard deviations of less than 1 mm to avoid spurious correlations with temperature and precipitation predictors (introduced in section 2b) and exclude snow on sea ice and permanent ice caps such as Greenland.
For comparison, we consider daily SWE, near-surface air temperature, and precipitation over 1979–2015 from NASA’s MERRA (Rienecker et al. 2011), based on version 5 of the Goddard Earth Observing System Model (GEOS-5) and version 5.2.0 of the Atmospheric Data Assimilation System (ADAS). MERRA provides a spatially complete gridded meteorological dataset through the assimilation of observations from balloons, aircraft, ships, buoys, and satellites. MERRA does not assimilate precipitation, except for microwave-retrieved rain-rate observations over ocean areas. These observations do not have a strong impact on MERRA’s precipitation, which is mostly influenced by three-dimensional humidity observations derived from the Special Sensor Microwave Imager (SSM/I) and the Advanced Microwave Sounding Unit-A (AMSU-A) on the Aqua satellite (Rienecker et al. 2011). Snow cover is represented in MERRA by the response to surface meteorological conditions of an intermediate-complexity snow scheme coupled with a catchment-based land surface model. This snow scheme has up to three snow layers describing snow accumulation, melting, refreezing, and other phenomenology related to the growth and ablation of the snowpack (Stieglitz et al. 2001).
MERRA is known to be among the most consistent reanalyses with observations (Liston and Hiemstra 2011; Lindsay et al. 2014). For example, among seven reanalyses products examined by Lindsay et al. (2014), MERRA’s monthly near-surface air temperature (1979–2009) ranked as having the smallest bias and highest anomaly correlation with observations in the high latitudes , whereas monthly precipitation (1989–2009) had the second-smallest bias and anomaly correlations of about 0.7 for every month. Also, Rapaić et al. (2015) showed that MERRA is among the two reanalyses with the closest overall agreement to the Canadian Gridded (CANGRID) Arctic-averaged annual temperature and precipitation series in 1980–2010. Sospedra-Alfonso et al. (2016) performed a detailed comparison of MERRA’s daily SWE with in situ observations at eight sites spanning a variety of climatic and land surface conditions in Canada. They found that MERRA represents the climatological seasonal cycle and interannual variability of SWE at these locations reasonably well, with correlation for annual maximum SWE averaged over the eight sites approximately equal to 0.6. Overall agreement with the in situ data was found comparable to, and in some cases better than, GlobSnow, a gridded SWE analysis that combines in situ data with SWE estimations from passive satellite microwave measurements (Takala et al. 2011) but is not available over mountainous regions. MERRA’s weaknesses include an anomalous increase in global mean precipitation after 1998 (Zhang et al. 2012), a positive cold season temperature and humidity bias over the Arctic (Serreze et al. 2012), a warm temperature bias over Eurasia during October (Brown and Derksen 2013), and a warm surface temperature bias during spring and summer (Jakobson et al. 2012). We consider MERRA data at their original ⅔° longitude × ½° latitude resolution.
b. Snow and meteorological indices
We employ snow and meteorological indices to characterize the interannual variability of SWE in a 30-yr window moving at 1-yr increments in 1850–2100. The snow indices, defined in a similar manner to Sospedra-Alfonso et al. (2015), include SWE monthly means in January (SWEJan), February (SWEFeb), March (SWEMar), and April (SWEApr) that represent the water content of the snowpack well after the snow onset and up to the start of the snowmelt in most of the mid- to high latitudes.
The meteorological indices are employed as predictors of the snow indices and are defined as follows. For every snow index, because snow amount is a function of meteorological conditions over some prior period, we define multiple predictors characterizing the aggregation of precipitation and temperature in previous months, starting at the beginning of the snow year. For example, SWEJan has six predictors of temperature and six predictors of precipitation given by , , , , , , , , , , , and denoting the average (total) temperature (precipitation) from August to January, September to January, and so on. Similarly, SWEFeb has seven predictors both for temperature and precipitation, and so forth.
c. Selection of best climatic predictors through partial correlation
For each grid cell included in our analysis, ensemble member, and sliding 30-yr interval in 1850–2100, we produce linearly detrended time series of annual values of each snow and meteorological index. We detrend the time series to distinguish correlations arising from interannual variability from those resulting from long-term trends, particularly in the future emission scenarios. We then use partial correlation to establish whether there is a linear relationship between monthly SWE and a given climatic predictor while controlling the influence of other predictors. In this way, we avoid “inflating” correlation values and attributing a relationship between SWE and one predictor that may be an indirect consequence of the others. For example, the partial correlation of SWEJan and gives a measure of the linear association between these indices excluding any variability that is linearly associated with that of all T predictors of SWEJan. This is done for every P predictor of a given snow index, and conversely partial correlations are computed between this snow index and every T predictor excluding variability that is linearly associated with all P predictors. We then select the meteorological indices (one for P and one for T) that have the highest partial correlation with the snow index as predictors for that snow index. These partial correlations for the best P and T predictors are denoted rP and rT, respectively. A high partial correlation indicates that the meteorological variable is strongly linearly associated with a particular snow index and thus can be considered a predictor of this snow index on physical grounds. If, for example, is the best P predictor of SWEJan for a given grid cell and 30-yr period, having a statistically significant partial correlation with this snow index that is higher than for the best T predictor, we conclude that precipitation is the primary driver of January mean SWE interannual variability at this location and time period. This is an indication that precipitation has been mainly in the form of snow and that temperature has remained sufficiently low to preserve the snowpack during October–January.
Note that because partial correlations exclude SWE variance attributable to the covariance of T and P that is present in the full correlations, effects of T predictors on SWE are quantified independently of effects of P predictors, and vice versa. Thus, for a given seasonal precipitation amount, the partial correlation with T will characterize how much of the precipitation falls as rain versus snow, since seasonal temperature will tend to determine the partition of precipitation into rain or snow. Conversely, for a given seasonal temperature, the partial correlation with P will characterize how much total precipitation, falling as rain or snow, influences SWE.
d. Skill of temperature and precipitation predictors of SWE variability
To further characterize temperature and precipitation as drivers of SWE variability, we estimate every snow index using multiple linear regression:
with fitting parameters α, β, and γ, where M denotes January, February, March, or April, and and are the best P and T predictors of as defined in section 2b, with and indicating the earliest aggregated month for P and T predictors, respectively, which are generally different.
To assess the statistical model in Eq. (1), we compute the anomaly correlation coefficient (ACC) and mean square skill score (MSSS) using determined through Eq. (1) to predict CanESM2 monthly values of . ACC is sensitive to the phase of estimated anomalies whereas MSSS penalizes phase errors, amplitude errors, and mean biases (Murphy 1988). Denoting as and the statistical and CanESM2 snow indices, respectively, for snow year and obtained from ensemble member , ACC in a Y-year window is computed as
where the primes denote anomalies with respect to the climatological values and σ denotes standard deviations. MSSS is computed using SWE climatology as a reference:
ACC and MSSS are computed from ensemble members in a years running window in 1850–2100.
e. Climatological seasonal temperature as a proxy for geographic location
To establish the geographical dependence of the relative influence of temperature and precipitation on SWE variability in the Northern Hemisphere, we employ the climatological seasonal temperature Ts of the aggregated months for the best P predictors as a proxy for geographic location. Specifically, for every snow-covered grid cell and snow index, we first compute for each year the average temperature over the aggregated months (i.e., season) of the best P predictor of this snow index (as defined in section 2c) and then the climatological value over each 30-yr moving window in 1850–2100. In this way, we assign a value of Ts to each snow-covered grid cell for each snow index and 30-yr period. We employ the terminology climatological “seasonal” temperature because the aggregated months of the best P predictors of a given snow index typically span, as we shall see in section 3b, from an early stage of the snow season up to the month specified by this snow index. The “season” lengths are thus variable in this context and generally are dependent on latitude/longitude and elevation.
f. Dependences on Ts and identification of a Ts threshold
For each snow index, considering all applicable grid cells for a given 30-yr window, we can characterize the partial correlations with the best P or T predictors as functions of climatological seasonal temperature Ts. As we show in section 3c, these functional relationships, which are generally sigmoidal with opposite dependences on Ts, characterize the geographic dependence of the relative effects of temperature and precipitation on snowpack variability. To determine the dependence of P partial correlations on Ts, we proceed as follows. For a 30-yr window, we first produce scatterplots of rP versus Ts from all ensemble members and applicable grid cells. We then partition the Ts axis into nonoverlapping bins (we consider bin sizes of 0.5°C or less) and compute the median of rP in every bin. The midpoint of every Ts bin is thus associated with the median of rP within that bin. The resulting values taken over all bins provide our estimate of the functional dependence of P partial correlations on Ts. We proceed similarly with the best T predictors, leading to two profiles of partial correlations, one for P predictors and one for T predictors, for each snow index. These profiles vary monotonically with opposing dependences on Ts and thus have a single intersection that defines a temperature threshold Tth separating two regimes: one where Ts > Tth and median rT > rP in which snowpack is mainly driven by temperature (T driven) and another where Ts < Tth and median rP > rT in which snowpack is mainly precipitation driven (P driven). To better quantify the likelihood of a P- or T-driven snowpack at a given location (i.e., for a given Ts), we also estimate the Ts interval corresponding to the intersection of rP and rT 33rd–66th-percentile ranges.
3. Results and discussion
a. Skill of the statistical model for SWE indices
We first assess the ability of the best T and P climatic predictors defined in section 2d to explain snowpack variability by means of the statistical model in Eq. (1). In doing so, we consider the ACC and MSSS skills defined in Eqs. (2) and (3) for each snow index and for the time periods 1850–79, 1981–2010, and 2070–99 in CanESM2 historical and RCP 8.5 scenario runs.
Figure 1 shows the geographic distributions of ACC. Generally, ACC is high , and tends to be highest in the continental interiors of Asia and North America and lower in coastal and peripheral snow-covered regions. ACC values are similar for each snow index, although somewhat lower in April as a result of the snowmelt and consequent decrease in snow-cover extent in the midlatitudes. ACC values for 1850–79 and 1981–2010 are very similar, although there are some reductions in the later period over the Tibetan Plateau and southern Rocky Mountains. By 2070–99, there are significant reductions in the area of appreciable ACC, mostly as a result of the northward retreat of the snow cover. However, ACC remains high in most of the snow-covered regions.
The geographic distribution of MSSS for each snow index and time period is shown in Fig. 2. MSSS is positive everywhere indicating that the error variance of Eq. (1) is always smaller than the variance of the predictand. As for ACC, MSSS is generally high and tends to be highest in continental interiors, with lower values in peripheral snow-covered regions. The values of MSSS are similar in 1850–79 and 1981–2010, although slightly reduced in the later time period. MSSS tends to be smaller in 2070–99, both in magnitude and extent.
The generally high values of ACC and MSSS indicate that the simple statistical model in Eq. (1) provides reasonably good estimates of SWE interannual variability for suitable T and P predictors regardless of emission scenario. Predicted SWE anomalies have generally the same sign as those of CanESM2 monthly SWE and tend to have the same magnitude, particularly in the higher latitudes and continental interiors. This shows that the combined effect of P and T predictors can be employed to characterize SWE interannual variability. In the following, we study the behavior of these T and P predictors and their influence on SWE variability separately to determine the regimes in which one type of predictor is less, equally, or more predominant than the other at driving SWE.
b. Behavior of the best T and P predictors
We examine here the geographic distribution and relation to climatological seasonal temperature Ts of the best T and P predictors for every snow index. The maps and box plots in Figs. 3–6 are obtained by computing median values of the number of aggregated months and climatological seasonal temperature Ts based on the five ensemble members at each grid cell. In the maps, we show the number of months between snow onset and the earliest aggregated month for best P predictors and separately the number of aggregated months of best T predictors. Snow onset is defined as the first day of at least 7 consecutive days with SWE > 4 mm. We exclude grid cells with negative median partial correlations between snow indices and best P predictors as these either have too little snow to establish robust statistical relationships or lack snow entirely. Analogously, we exclude grid cells with positive ensemble median partial correlations between snow indices and best T predictors, which can occur either if there is little snow or when mean temperatures are sufficiently cold such that T predictors and monthly SWE variability are insufficiently related during the snow accumulation period. The number of aggregated months in Figs. 4–6 are counted from the month associated with the snow index toward earlier months.
Figure 3 shows the geographic distribution of the number of months between snow onset and the earliest aggregated month of best P predictors in 1981–2010 and 2070–99 from CanESM2 historical and RCP 8.5 scenario runs. Overall, P predictors tend to have aggregations going back to the earliest stages of the snow season regardless of snow index, with earliest months typically coincident with, or one month later than, the timing of snow onset. In 1981–2010, snow onset occurs in October or November in most of the mid- to high latitudes (not shown). Aggregations starting two or more months after snow onset are rare, although they occur in some portions of Europe, the coastal regions, and/or portions of the snow-cover periphery where there may be significant snowmelt events during the snow accumulation period. In 2070–99, aggregations start typically about one month after snow onset, which is in November or December in most of the mid- to high latitudes (not shown). The shift toward later aggregations relative to snow onset is most likely the result of a delayed snow season combined with the coarse temporal resolution employed to determine the aggregation of the predictors. This may also explain, at least in part, the slight increase of regions with aggregations starting two or more months after snow onset during this period.
The geographic distribution of the number of aggregated months of best T predictors is shown in Fig. 4. Overall, T predictors tend to have short aggregations starting in the months closest to those of the snow indices. In 1981–2010, a large portion of the midlatitudes has best T predictors for SWEJan and SWEFeb starting at or later than November. For SWEMar and SWEApr there is a tendency for aggregations of T predictors to start in February or March, or mostly April in the latter case. These tendencies are largely unchanged in 2070–99, except that the extent of the snow-covered regions is generally reduced. Snow-covered regions with T predictors having positive correlations with snow indices are not colored.
Figure 5 shows box plots of Ts versus the number of aggregated months of the best P predictors of every snow index in 1981–2010 and 2070–99. The values of Ts tend to decrease with increasing number of aggregated months, until aggregations through October (1981–2010) or November (2070–99) regardless of snow index. Aggregations up to October–November (1981–2010) and November–December (2070–99) encompass the great majority of grid cells (typically %). This indicates that snowpack “memory” of previous snowfalls tends to increase with colder mean seasonal temperature and goes back to the early stages of the snow season, consistent with the maps in Fig. 3. The few locations with P predictors having more aggregated months are either in the higher latitudes where snow seasons start in August or September or in regions of ephemeral snow where Ts is typically high and correlations are not statistically significant. Two important changes in 2070–99 relative to 1981–2010 evident in Fig. 5 that are a consequence of a warming climate are 1) the marked increase in the values of Ts percentiles (5°–10°C warmer for the minima, over the number of aggregated months, of median Ts) and 2) a shift toward fewer aggregated months in the best P predictors resulting from later snow onsets. This indicates that by the end of the century, SWE variability will likely be driven by warmer and seasonally shorter P predictors. This shift toward seasonally shorter P predictors is already evident in 1981–2010 when compared to similar plots for early periods, particularly 1850–79 (not shown), as the largest percentage of grid cells in 1850–79 have P predictors with aggregations up to October, as opposed to November in 1981–2010.
Box plots for the number of aggregated months of best T predictors of every snow index in 1981–2010 and 2070–99 are shown in Fig. 6. In this case, the dependence of Ts and the number of aggregated months is less systematic than in Fig. 5, as Ts is based on the aggregated months of the best P predictors. For SWEJan and SWEFeb, the level of aggregations is more or less equally distributed among grid cells with a tendency for short aggregations. For SWEMar and SWEApr, short aggregations predominate, suggesting that T predictors of SWE variability are influenced mostly by recent snowfalls and/or snowmelt events during the late stages of the snow accumulation period and/or start of the snowmelt season, which is consistent with Fig. 4 and in sharp contrast with the long memory of P predictors shown in Figs. 3 and 5.
c. Dependences on Ts and determination of a Ts threshold
As expected from the properties of the climatic predictors discussed in the previous section, the effects of temperature and precipitation on snowpack variability should depend on climatological seasonal temperature and therefore on geographic location. We examine the dependences of the strength of the linear relationship between snow indices and the best P or T predictors with Ts and then show that there is a Ts isotherm, largely independent of snow accumulation period and emission scenario, that divides the snow cover into two regions: one where the snowpack is predominantly driven by precipitation and the other where it is mainly driven by temperature.
The comparative influences of P and T on snowpack variability as a function of Ts formulated as described in section 2f are shown in Fig. 7. Here, each panel consists of scatterplots of Ts versus the partial correlations of each snow index with its best P or T predictors in 1850–79 and 2070–99 from historical and RCP 8.5 scenario runs, emphasizing the behavior of a central tercile range between the 33rd and 66th percentiles. Overall, the relationships between rP or rT and Ts are nonlinear over the range of climatological temperatures. In particular, rP values tend to increase as Ts decreases from approximately 0° to −10°C and remain roughly constant for colder and warmer temperatures. By contrast, rT absolute values tend to decrease with decreasing Ts from approximately 0° to −20°C and remain roughly constant outside this range, although in 2070–99 there are almost no regions with Ts < −20°C as a result of projected warming. Tracking the median rP and rT values with respect to Ts allows us to estimate a threshold temperature Tth representing a crossover point below which the snow season is sufficiently cold such that the snowpack is mostly driven by precipitation. Above that threshold, temperature tends to drive the snowpack. In each panel of Fig. 7 Tth lies between −10° and 0°C.
2) Stability with respect to time and emission scenario
Figure 8 shows for every snow index the time evolution of Tth and the temperature band derived from the intersection of the central T and P partial correlation terciles based on successive 30-yr windows in 1850–2100 from historical and RCP 2.6, 4.5, and 8.5 scenarios. For each snow index, the temperature threshold is largely independent of time and emission scenario, with little linear trend (mostly <0.01°C per period of given emission scenario) and standard deviation typically ≤0.6°C. There is a modest progressive decrease of Tth during the snow season, from mean values ranging from −4.1°C for SWEJan to −6.7°C for SWEApr, differing little between the historical period and future scenarios.
For each snow index, the relative stability of Tth with time and emission scenario results from the near time invariance of the dependence of rP and rT on Ts. Figures 9 and 10 show, respectively, the differences of the median rP and rT profiles during 1850–2100 relative to 1850–79. The statistically significant differences at p < 0.01 according to the Mood’s median test are shown in color, and they represent 3.3% (2.7%), 3.7% (5.5%), 3.4% (2.4%) and 5.4% (3.4%) of all differences shown for P (T) profiles corresponding to SWEJan, SWEFeb, SWEMar, and SWEApr, respectively. Overall, these differences are relatively small for most values of Ts and tend to be slightly negative for future times, mostly for P profiles, implying a slight decrease in the partial correlations. For rT, differences are most notable between −20° and −10°C, particularly for SWEJan and SWEFeb. The largely invariant profiles of rP and rT with time are at least partly a consequence of employing the best P and T predictors for the snow indices, as these predictors have variable levels of aggregation depending on emission scenario. If the predictors’ level of aggregations were held fixed for different emission scenarios, then rP and rT median profiles would change more over time (not shown). This suggests that the invariance of P and T correlation profiles might be strengthened if the levels of aggregation in the P and T predictors were refined (e.g., counting aggregations in weeks or days instead of months).
3) Dependences of T- or P-driven snowpack variability on Ts
Because the partial correlation profiles and thus Tth remain relatively unchanged regardless of emission scenario, we can combine data from all grid cells, ensemble members and 30-yr moving windows in 1850–2100 in a single scatterplot to obtain a more accurate view of the relative influence of T and P on snowpack variability as a function of Ts, and thus a more accurate estimate of Tth.
Following this approach, we show in Fig. 11 the 33rd–66th percentiles, with the medians highlighted, of the partial correlations of each snow index with its best P and T predictors from CanESM2 historical and RCP 8.5 scenario runs in 0.25°C Ts bins. For each snow index, the median rP values are typically high and nearly constant for Ts < −12°C (or a slightly lower temperature for SWEApr). These values decrease approximately linearly to about 0.4 for Ts between −10° and ±2°C, and gradually decrease further to around 0.3 for higher Ts. On the other hand, the median absolute values of rT tend to remain low for Ts < −20°C, increase approximately linearly to about 0.8 at Ts±2°C, and remain almost constant for higher Ts. The relationship between the median values of rP or rT with Ts resembles a sigmoid, which can be fitted with a logistic function of the form
where a, b, c, and d are fitting parameters given in Tables S1 and S2. In particular, a and b are the lower and upper asymptotes of the sigmoid, respectively, c gives the steepness of the curve between asymptotes, and d is the Ts locus of maximum slope. For each fitted curve, the coefficient of determination is above 0.98, implying an excellent agreement with each of the rP and rT median profiles. As shown in the tables, the fitting parameters are largely invariant with respect to emission scenario.
The strength of the linear relationship between P predictors and snow indices, as well as its dependence on Ts shown in Fig. 11, is largely explained by the fraction of days with mean temperature below freezing in the aggregated months of the best P predictors. This is shown in Fig. 12 by pooling the five ensemble members in a single plot and, for each nonoverlapping bin in Ts of a given size (0.5°C in the figure), computing the median number of days with mean temperature below 0°C in the aggregated months of the best P predictors for each 30-yr window in CanESM2 historical and RCP 8.5 scenario runs. Overall, the fraction of days within a season having mean temperatures below freezing decreases with increasing Ts, in a similar manner to rP in Fig. 11. For Ts below −12°C (slightly lower for SWEApr) the fraction of days with temperature below freezing is typically large (>95%), which explains the saturation of high rP values for this climatological temperature range (Fig. 11). For Ts above −12°C, the fraction of cold days decreases with increasing Ts following a dependence that largely explains that in Fig. 11 and tends to saturate for Ts values well above freezing. As implied by Tth values, the regions with P-driven snowpack tend to have mean daily temperature below freezing at least 75% of the snow season, whereas T-driven snowpack tends to occur in regions where temperature is above freezing at least 25% of the snow season.
4) Climatological seasonal temperature threshold
The intersection of the P and T median correlation profiles in Fig. 11 defines Tth, which is given in Table 1 with its 33rd–66th-percentile range (defined at the end of section 2) as well as the value from the intersection of the analytical fits deduced from Eq. (4), for several sizes of temperature bins and emission scenarios. Overall, Tth changes little with emission scenario and shifts approximately 2°C toward colder temperatures between January and April, with the biggest change taking place from March to April (over 1°C). This is likely the result of increased monthly mean temperature from January to April, particularly from March to April, which coincides with the beginning of the snowmelt in large areas of the midlatitudes.
Figure 13 shows maps of Ts and Tth, the latter derived with the analytical fits in Eq. (4), for SWEJan, SWEFeb, SWEMar, and SWEApr in 1850–79, 1981–2010, and 2070–99 from the historical and RCP 8.5 scenario. In regions with Ts < Tth, SWE variability is mostly P driven, whereas in regions with Ts > Tth, SWE variability is mostly T driven. The Ts values differ relatively modestly between 1850–79 and 1981–2010, although notable changes occur in the Russian Arctic, where regions with Ts < −20°C have largely receded, and in the Tibetan Plateau and Rocky Mountains, where mostly isolated regions of P-driven SWE variability delimited by Tth are also reduced. As the snow-climate regions with Ts < Tth recede, the extent of the colder portion of the snow cover that is mostly driven by precipitation diminishes, much as the receding 0°C isotherm reduces the snow-cover extent. As a result, increasingly more northern regions and elevation zones enter the Ts > Tth regime, becoming more vulnerable to warming temperatures. In 2070–99, most of these regions have seasonal mean temperature above freezing as a result of the significant northward shift of the isotherms, leading to a widespread reduction in the snowpack and snow-cover extent, most notably in North America, Europe, and southern Asia. In this time period, P-driven snow-covered regions have become concentrated mostly northward of 60°N and, before April, in reduced ranges southward of 60°N in eastern Russia, Mongolia, northern Canada and China, and the Tibetan Plateau. In April, most snow climate regions southward of 60°N have become T driven, including most of the Tibetan Plateau.
5) Comparison with reanalysis
Comparisons of P and T partial correlation profiles for CanESM2 relative to the MERRA reanalysis in 1979–2015 are shown in Fig. 11. For SWEJan, SWEFeb, and SWEMar there is generally good agreement between CanESM2 and MERRA, with differences mostly for T-correlation profiles. For P-correlation profiles, this agreement is explained, at least in part, by the similar behavior between CanESM2 and MERRA of the fraction of days having mean temperature below freezing in the aggregated months of best P predictors (Fig. 12), as a higher fraction of sufficiently cold days during the snow season leads to a stronger association between SWE and precipitation variability. This agreement occurs despite differences in the aggregation levels of best P predictors in CanESM2 and MERRA. For example, maps of the earliest aggregated month of these P predictors in 1981–2010 (not shown) present similar patterns for both data products, but MERRA P predictors tend to have more aggregated months than CanESM2 in the higher latitudes (>60°N) and fewer aggregated months in regions of peripheral snow and the Tibetan Plateau. Overall, MERRA has stronger linear association between snow indices and best T predictors than CanESM2 for seasonal temperature below Tth, particularly for Ts between −20° and −10°C (Fig. 11). This is consistent with the higher rates of snowmelt in MERRA relative to historical simulations of CanCM4, which shares physical components with CanESM2, during August–March in 1981–2010 (Fig. 8 in Sospedra-Alfonso et al. 2016), as increased snowmelt should be an indication of a stronger association between SWE and temperature variability.
For SWEApr, MERRA typically has lower P correlations than CanESM2 for Ts between −20°C and Tth and higher absolute T correlations, which indicates that temperature variability plays a more significant role at driving mean April SWE in MERRA than in CanESM2 for such relatively cold mean temperatures. This is consistent with the early melt of Arctic snow in MERRA relative to historical simulations of CanCM4 (Fig. 5 in Sospedra-Alfonso et al. 2016), as well as the generally higher rate of snowmelt mentioned above, despite MERRA having colder mean temperature over snow than CanCM4 (Fig. 9 in Sospedra-Alfonso et al. 2016). A factor that may influence, at least in part, the reduced P correlations in MERRA for Ts < Tth is its significantly higher rate of sublimation relative to CanCM4 during January–April (peaking in March–April) (Fig. 8 in Sospedra-Alfonso et al. 2016). Snow sublimation may degrade the linear association between SWE and snowfall as it can lead to significant snow removal in dry cold winters (e.g., Liston and Hiemstra 2011). Another potential source of degradation in the partial correlation of snow indices and P predictors is snow loss by transport (e.g., Pomeroy and Gray 1994), including sublimation of blowing snow, but neither MERRA’s land surface model nor the CLASS 2.7 snow scheme account for snow redistribution by wind.
d. Case studies in North America
Because Ts may increase considerably in a given location as a result of a warming climate, the local strength of the linear relationship between snow index and P or T predictors may change significantly according to the sigmoidal curves in Fig. 11. As implied by Figs. 11 and 12, a snow accumulation period having over 85% of days with mean temperature below 0°C (as inferred from the lowest end in the temperature threshold interval derived with the partial correlations' 33rd–66th-percentile range) is likely to produce a P-driven snowpack. Conversely, the degree to which daily mean temperature 85th percentiles during the snow accumulation period exceed freezing should provide an indication of the degree to which temperature drives snowpack variability. To examine this, we show in Fig. 14 the time evolution of the ensemble median of Ts and , as well as the partial correlations of SWEApr with its best P and T predictors, from CanESM2 historical and RCP 8.5 scenario runs, for four locations in western North America. We consider locations with initially different values of Ts < 0°C and above or below freezing.
The location depicted in Fig. 14a is at the lowest latitude and has the highest seasonal mean temperature among all considered. At all times, this location has Ts above Tth and well above 0°C, resulting in a T-driven snowpack with moderate to low values of rP and higher absolute values of rT . This location becomes snow-free in April by 2030 when Ts exceeds the freezing level and rP vanishes. For the other locations, Ts is initially below Tth whereas is initially below freezing, and they remain so until the end of the first half (Fig. 14b) and the second half (Figs. 14c,d) of the twenty-first century. In these time periods, rP is largely stable about 0.8–0.9 and starts declining when consistently approaches 0°C. Conversely, absolute rT is moderate to low during these initial time periods and tends to increase as increases, leading to a crossover from P- to a T-driven snowpack when Ts reaches the temperature threshold Tth and reaches the 0°C isotherm.
4. Summary and conclusions
We have studied the relative effects of temperature and precipitation on snowpack variability in the Northern Hemisphere in CanESM2 historical (1850–2005) and future (2006–10) scenario climate simulations. Snowpack variability, characterized here by monthly SWE values in January, February, March, and April, is typically explained in terms of predictors consisting of accumulated precipitation from the early stages of the snow season and/or the average temperature of the months closest to that of the snow indices. Consequently, there is no single pair of best P and T predictors of snowpack variability for the whole Northern Hemisphere but rather differing predictors with different levels of aggregation depending on geographic location, time period, and emission scenario.
Even though the response of snowpack to T and P variability at a given location is generally strongly dependent on emission scenario (Fig. 14), the overall response of the snowpack to predictors of a given climatological seasonal temperature was found to be largely independent of time and emission scenario (Figs. 7, 9, and 10). Thus, we found that there is an isotherm mostly characterizing the thermal properties of the snow accumulation period (−5±1°C depending on snow index and largely independent on emission scenario), below which precipitation primarily drives snowpack variability and above which temperature is the main driver. Moreover, the strength of the linear relationship between monthly SWE and best P or T predictors depends on the climatological temperature of the snow accumulation period Ts according to a sigmoidal function.
We found three distinct thermal regimes of the snowpack response to temperature variability determined by this sigmoidal function: (i) Ts < −20°C, where SWE is nearly insensitive to temperature variability; (ii) −20°C < Ts < ±2°C, where SWE becomes increasingly sensitive to temperature variability (approximately linearly) with increasing Ts; and (iii) Ts > ±2°C, where SWE is highly sensitive to temperature variability, nearly independently of Ts. Likewise, we found three distinct thermal regimes of snowpack response to precipitation variability: (i) Ts < −12°C, where SWE is highly sensitive to precipitation variability nearly independently of seasonal temperature; (ii) in the range of −12°C < Ts < ±2°C, where SWE sensitivity to precipitation variability tends to decrease approximately linearly with increasing Ts; and (iii) Ts > ±2°C, where SWE sensitivity to precipitation variability is weak and nearly independent of Ts. The transition from strong to weak sensitivity of SWE to precipitation variability is largely explained by the fraction of days with temperature below freezing during the snow accumulation period (Fig. 12). We found that regions having typically more than 25% of the snow accumulation period with daily mean temperature above freezing tend to have snowpack amount largely driven by temperature and are thus more vulnerable to a warming climate regardless of emission scenario.
Using this methodology we can assess snowpack sensitivity to variations in precipitation and temperature, bearing on whether snow amounts are likely to increase, decrease, or stay about the same as a result of a changing climate. In a warming climate, isotherms shift northward diminishing the precipitation-driven portion of the snow cover, much as the snow cover itself diminishes with the shift of the 0°C isotherm. Increased susceptibility of SWE to temperature variability is thus expected in snow-climate regions with Ts about or above −5°C, which are mostly located in the coastal and more southerly portion of the Northern Hemisphere midlatitudes in the present-day climate. By contrast, SWE in regions with Ts well below −20°C mostly located in northeastern Siberia and the Canadian Archipelago in present-day climate is largely unaffected by temperature variability, implying that increasing precipitation could lead to increasing snow accumulation in those regions. These results are consistent with findings by Räisänen (2008) and Brown and Mote (2009).
We showed that the best P predictors of SWE variability have time aggregations ranging typically from the start of the snow accumulation period up to the month considered (Figs. 3 and 5), whereas the best T predictors tend to be determined by the previous 1 to 2 months from the month considered (Figs. 4 and 6), particularly for March and April, which encompass peak SWE for most of the mid- to high latitudes. This rule of thumb can serve as the basis of a simpler implementation of the methodology presented here, as it avoids the necessity of considering all P and T predictors of SWE at a given location, although optimal results are obtained by using P and T aggregations based on local climatologies as done here. A potential application of this simpler approach is estimating SWE anomalies based on such P and T predictors. Indeed, preliminary results indicate (not shown) that the regression coefficients α and β in Eq. (1) for standardized predictors and predictand (i.e., P, T, and SWE having climatological means removed and normalized with the standard deviations) follow similar sigmoidal relationships with Ts as those in Fig. 11 for the partial correlations. Then, by using the analytic fit in Eq. (4) to such sigmoidal relationships, together with Ts as well as P and T predictors derived with the rule of thumb mentioned above, we can use Eq. (1) to estimate SWE anomalies based only on local temperature, precipitation, and climatology of the timing for snow onset. Such application is outside the scope of this work, and we leave it for a separate study.
The analysis presented here has provided insight into the snowpack response to climate variability and change in the Canadian Earth system model CanESM2. We have introduced a framework that can be employed to study the roles of temperature and precipitation in snowpack variability as a function of latitude/longitude and elevation, as well as time. In addition, the sigmoidal relationship between climatological seasonal temperature and the strength of the linear association of snow indices and P or T predictors (Fig. 11), as well as the relative time invariance of the temperature threshold separating the snow cover into P- or T-driven regions (Fig. 8), can be used as a model diagnostic that can serve as a tool to assess the impact of a warming climate on the snow cover. As hypothesized by Räisänen (2008), SWE increase in the coldest portions of Eurasia and North America combined with the widespread SWE decrease in milder regions suggests a strong dependence between projected SWE changes and present-day winter temperatures. We have assessed the functional dependence of this relationship (Fig. 11) and found that future SWE variability in any given snow-climate region are largely dependent on its present-day mean snow-seasonal temperature through the correlations of P and T predictors with snowpack variability determined by this function. This functional dependence could further be employed to improve SWE projections by correcting for biases in simulated present-day temperature, as suggested by Räisänen (2008).
We are grateful to Ross Brown for useful comments on an earlier draft of this paper and for pointing out relevant publications. We are also grateful to two anonymous reviewers and the editor for insightful comments that helped improve the presentation of this work. This research is partly supported by the Natural Sciences and Engineering Research Council of Canada’s (NSERC; Grant 476404) Climate Change and Atmospheric Research (CCAR) program (Grant 433874) through the Canadian Sea Ice and Snow Evolution (CanSISE) Network (www.cansise.ca). R. S.-A is funded by an NSERC Visiting Fellowship.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-16-0612.s1.