Change in snow phenology and its potential feedback to temperature in the Northern Hemisphere over the last three decades

Table 1 summarizes the datasets used in this study. Daily snow depth and surface air temperature data from August 1979 to July 2006 in China, Europe, Russia, Canada and USA were obtained from the National Meteorological Information Center of the China Meteorological Administration (Peng et al 2010), European Climate Assessment & Dataset (ECA&D) (available at http://ecad.knmi.nl/dailydata/predefinedseries.php), Russian Research Institute of Hydrometeorological Information—World Data Center (RIHMI-WDC) (Bulygina et al 2011), the National Climate Data and Information Archive of Environment Canada (available at http://climate.weatheroffice.gc.ca/prods_servs/index_e.html) and the United States Historical Climatology Network (USHCN) (available at http://cdiac.ornl.gov/epubs/ndp/ushcn/ushcn.html), respectively. The albedo dataset used in this study was from EOS-Terra-MODIS MOD43C3 Nadir BRDF-Adjusted Reflectance (NBAR) 2001–2006 products, at a spatial resolution of 0.05° × 0.05° and 8 d interval including black sky albedo (BSA) and white sky albedo (WSA) over short-wave broadband (0.3–5.0 μm) (Schaaf et al 2002). WSA shows similar results to BSA, so we used BSA to illustrate the results here. The bias in MODIS albedo data (Collection V005) is mostly less than 5% (Liu et al 2009). The incoming short-wave radiation dataset used in this study was from ERA-Interim reanalysis produced by European Centre for Medium-Range Weather Forecasts (ECMWF) over the period 2001–2006, at ∼0.7° ×∼ 0.7° spatial resolution and 6 h temporal resolution, with the bias of ∼2 W m⁻² (Dee et al 2011).

To avoid the impact of ephemeral snow on the ground on the snow cover duration, the onset date of snow cover for year (t) at each meteorological station is defined as the first continuous five days on which snow was observed to cover the ground surface in the fall of the year (t − 1); the end date of snow cover is defined as the last continuous five days when snow cover was observed in the spring of the year (t). The snow cover duration is defined as the number of days from the onset date of snow cover to the end date of snow cover for the year (t). For any site, the year (t) was omitted from the analysis if there were one or more days without snow depth data during the period from 1 August of year (t − 1) to 31 July of year (t). Furthermore, in order to reduce uncertainties caused by the short time period of observation, we consider only stations for which at least 20 years of daily snow depth data are available during the study period. Among 636 meteorological stations, which have at least 20 years of daily snow depth data, only 575 stations have at least 20 years of both daily snow depth and monthly air temperature data during the study period 1980–2006. Maximum snow depth during snow season is calculated as the maximum snow depth record from the onset date of snow cover to the end date of snow cover. We referred spring and fall as the snowmelt (February–July) and onset (September–January) season.

To map the spatial distribution of the trend of the end date of snow cover (D_E), we applied the trend detection methodology from Weatherhead et al (1998) which consider the autocorrelation of the underlying noise processes (assuming noise is autoregressive of the order of 1), at each of the 636 meteorological stations with at least 20 years of daily snow depth data during the period 1980–2006. We similarly calculated the trend of the onset date of snow cover (D_O), the duration of snow cover and the temperature at each site.

For the sensitivity of D_E to temperature, to remove the effects of snow depth on D_E, we regressed the end date D_E as dependent variable with average air temperature during the month before snowmelt (T_M−1) and maximum snow depth (SD_M) as independent variables.

$$\begin{eqnarray} &&{D}_{\mathrm{E}}=a\times {T}_{\mathrm{M}-1}+b\times S{D}_{\mathrm{M}}+c \end{eqnarray} \tag{ 1 }$$

where a and b are regression coefficients, and c is the residual error. D_E mainly depends on the temperature before snowmelt and the maximum snow depth to melt (see section 3). We defined the regression coefficient of T_M−1 (a) as the sensitivity of D_E to temperature.

To examine the effect of positive feedback of earlier snowmelt on temperature, we performed a multiple linear regression (equation (2)) using average air temperature during snowmelt month (T_M) as the dependent variable, and D_E and average air temperature during the month before snowmelt (T_M−1) as independent variables.

$$\begin{eqnarray} &&{T}_{\mathrm{M}}=s\times {D}_{\mathrm{E}}+p\times {T}_{\mathrm{M}-1}+q \end{eqnarray} \tag{ 2 }$$

where s and p are regression coefficients, and q is the residual error. We assumed that interannual variations of T_M come from general atmosphere circulation and the local snow albedo feedback, which are reflected by the term p × T_M−1 and s × D_E, respectively. The slope of the regression, s, defines the sensitivity of temperature to D_E.

Similarly, for the sensitivity of temperature to D_O, we performed the multiple linear regression as equation (2), but using the average air temperature during snow onset month (T_O) as the dependent variable, and D_O and average air temperature during the month before snow onset (T_O−1) as independent variables as equation (3)

$$\begin{eqnarray} &&{T}_{\mathrm{O}}=s\times {D}_{\mathrm{O}}+p\times {T}_{\mathrm{O}-1}+q, \end{eqnarray} \tag{ 3 }$$

where s and p are regression coefficients, and q is the residual error. The slope of the regression, s, defines the sensitivity of temperature to D_O.

For the incoming short-wave (SW) radiation at surface at D_E or D_O in each station, we used average daily incoming SW radiation at D_E or D_O during the period 2001–2006 in the pixel (6 h, ∼0.7° ×∼ 0.7°) where the station located. For the surface albedo, we resample the station surface albedo from MODIS MOD43C3 (8 d, 0.05° × 0.05°). For each station, we used maximum and minimum black sky albedo (BSA) and white sky albedo (WSA) from MODIS MOD43C3 during the month of snowmelt as the surface albedo before and after snowmelt (figure S1 available at stacks.iop.org/ERL/8/014008/mmedia), respectively. Then we took the difference in BSA before and after snowmelt as δBSA. The product of incoming SW radiation during snowmelt (snow onset) by δBSA is the difference in absorbed short-wave energy (ASR) between snow free and snow covered surface (δASR).

Figures 1(A) and (B) show the spatial patterns of mean onset and end date of snow cover during the period 1980–2006. There are clear latitudinal gradient patterns for the onset (earlier with higher latitude) and end (later with higher latitude) date of snow cover (figures 1(A) and (B)). Then the mean snow cover duration shows longer snow duration period with higher latitude (figure 1(C)). We first applied the trend detection method (Weatherhead et al 1998) to the time series of the onset and end date of the snow cover in each year, and of the snow cover duration at each station (figure 2). A decrease in snow cover duration is found for 63% of the stations (figure 2(C)). This result agrees with the observation from Scanning Multi-channel Microwave Radiometer (SMMR) satellite data with global coverage at 25 km resolution and 2 d revisit (Choi et al 2010, Takala et al 2009). The rate of shrinking of snow duration is greater than five days per decade at more than 34% of the 636 stations (figure 2(C)). Interestingly, North America and Eurasia show different trends for the end date of snow cover in spring. A trend towards an earlier end date of snow cover is found at 70% of the Eurasian sites, with 28% of the sites showing an advance of snowmelt in spring as fast as five days per decade. By contrast, the trends in the end date of snow cover over North America exhibit a fragmented spatial pattern. At the continental scale, the end date of snow cover over the last 27 years shows a rapid advance for Eurasia (−2.6 ± 5.6 d decade⁻¹) and virtually no trend for North America (0.1 ± 5.8 d decade⁻¹). The onset date of the first snow shows a trend towards later dates in Eurasia (1.3 ± 4.9 d decade⁻¹) and North America (1.1 ± 4.9 d decade⁻¹). These results are consistent with the trends of vegetation phenology estimated by satellite observation and some terrestrial biogeochemical cycle models (Smith et al 2004, Tucker et al 2001, Piao et al 2008, Jeong et al 2011). The green-up of vegetation in spring for many ecosystems north of 25°N was modeled to follow a trend towards earlier dates in Eurasia (Piao et al 2008, Jeong et al 2011), in agreement with satellite observation of vegetation reflectance (Zhou et al 2001).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We hypothesize that the geographically different trends in snow cover phenology between North America and Eurasia for the end date of snow cover, are associated with local temperature trends. From co-located temperature records at 575 snow measurement stations (see datasets and methods), we found that in the seasons of snow melt (February–July), 85% of the Eurasian meteorological stations indeed show an increase in temperature during the month of snowmelt (T_M) and T_M increased by 0.21 ± 0.18 °C decade⁻¹ for all Eurasian stations, whereas only 49% of the North American stations show an increase in T_M and at a smaller rate of 0.00 ± 0.21 °C decade⁻¹ for all North American stations. During the fall (seasons of snow accumulation beginning), however, air temperature appears to increase at a similar rate over North America (0.10 ± 0.15 °C decade⁻¹) as over Eurasia (0.13 ± 0.14 °C decade⁻¹), in agreement with the hypothesis of a co-variation between temperature and the presence of snow.

In addition to spring air temperature, we also hypothesize that the end date of snowmelt is also related to snow depth, a proxy for the amount of snow to be melted (Stone et al 2002). Among the 636 stations in this study, 95% of them show a positive correlation of the end date of snowmelt with the maximum snow depth during the snow season (statistically significant at 49% of the 636 stations). This corroborates the work by Stone et al (2002) who reported that snow water equivalent influences the snowmelt date at Barrow station. We found that the spatial distribution of the trend in maximum snow depth during the snow season can only explain a limited (7.0%), but statistically significant (P < 0.001) portion of the observed distribution in the trend in end date of snowmelt over the 636 stations.

To estimate the coupling between spring snow disappearance and air temperature, we first performed a correlation analysis between the end dates of snow cover each year and T_M−1 (monthly air temperature during the month before snowmelt) across 575 stations where at least 20 years of daily snow cover and collocated air temperature data were available during the study period. A negative correlation between the end date of snow cover and T_M−1 was obtained at more than 94% of the 575 stations (statistically significant at 40% of the stations), suggesting that spring air temperature plays a critical role in controlling the timing of snowmelt. Because both maximum snow depth during snow season and spring air temperature can influence the end date of snow cover, equation (1) was applied to estimate the sensitivity of the end date of snow cover to T_M−1 after controlling the effect of snow depth on the end date of snow cover for each station (partial correlation analysis). We found that, for a given increase of T_M−1 of 1 °C, the end date of snow cover advances on average by 1.6 ± 1.8 d (figure 3(A)). In addition to the impact of warmer spring temperatures on the end date of snow cover, a feedback of snowmelt on air temperature is also expected (Groisman et al 1994, Déry and Brown 2007, Hall and Qu 2006). The magnitude of this feedback is analyzed below.

Zoom In Zoom Out Reset image size

In this section, we made an attempt to quantify the feedback of snowmelt on spring temperature. To do so, we established a relationship between the end date of snow cover (D_E) and temperature during the month of snowmelt (T_M) for each station (retaining only sites with at least 20 years of daily snow cover and air temperature data from 1980 to 2006). We found that 98% of the stations exhibited a negative correlation between T_M and D_E (statistically significant at 67% of stations). At 69% of the stations, T_M shows an even stronger negative correlation with D_E than with T_M−1. We performed multiple linear regression analysis using T_M as the dependent variable and both D_E and T_M−1 as independent variables (predictors) to further remove the effect of autocorrelation of temperature between two consecutive spring months (see datasets and methods). The result was that changes in D_E are significantly and negatively correlated with T_M at 332 of the 575 NH stations, when controlling for the effect of T_M−1. The probability that this occurs purely by chance is estimated to be less than 0.1% if all station records are assumed to be independent. The linear regression slope, or sensitivity, of spring temperature to D_E is −0.077 ± 0.067 °C d⁻¹ (figure 3(B)). This sensitivity can be interpreted as a measure of the snow-temperature feedback. In other words, an earlier disappearance of snow cover by 1 d could result in an additional warming of monthly temperature of 0.077 ± 0.067 °C.

Likewise, in the fall (onset of snow) we found positive correlations between the onset date of snow cover (D_O) and the temperature during the month of snow onset (T_O) in 98% of the stations (statistically significant at 54% of the stations). In other words, a colder autumn coincides with an earlier snow onset. When we remove the effect of temperature during the month preceding snow onset (T_O−1), the sensitivity of T_O to D_O is 0.068 ± 0.051 °C d⁻¹ (figure 3(C)), a little smaller in absolute value than the sensitivity of T_M to D_E (−0.077 ± 0.067 °C d⁻¹). This result suggests a small feedback of snow cover on local temperature in autumn, a period of the year when the incoming solar radiation (the total amount of energy whose availability is controlled by albedo) is smaller than that in snowmelt season (spring).

In addition to temporal relationships at each site, we also found a negative (positive) correlation between trends in the end (onset) date of snow cover and trends in spring (autumn) temperature across the sites. We regressed the trends in T_M (T_O) against the trends in D_E (D_O) across the 575 stations, respectively (figure 4). The key result is that the trends in T_M (T_O) are negatively (positively) associated with the trends in D_E (D_O) across the sites with R =− 0.27, P < 0.001 (R = 0.22,P < 0.001). After removing the effect of the trend in T_M−1 (T_O−1), there is still a significant residual negative partial correlation between the trend in T_M (T_O) and the trend in D_E (D_O) across the stations (R = 0.20,P < 0.001;R =− 0.20,P < 0.001; figure 4). The spatial sensitivity of temperature to D_E (D_O), given by the linear slope between the T_M (T_O) and D_E (D_O) trend is −0.007 °C d⁻¹ (0.008 °C d⁻¹), which is smaller than the temporal sensitivity established above (−0.077 ± 0.067 °C d⁻¹ for the sensitivity of T_M to D_E and 0.068 ± 0.051 °C d⁻¹ for the sensitivity of T_O to D_O, respectively).

Zoom In Zoom Out Reset image size

The sensitivity of T_M (T_O) to D_E (D_O) depends among others, on the amount of extra absorbed incoming short-wave (SW) radiation due to earlier snowmelt (later snow onset), which is the product of incoming SW radiation during snowmelt (snow onset) and the difference of albedo between snow free and snow covered surfaces. We quantified two important driving factors that are hypothesized to control the sensitivity of T_M (T_O) to D_E (D_O): short-wave incoming radiation, and surface albedo. Firstly, we correlated the sensitivity of T_M to D_E to daily short-wave incoming radiation across the 575 selected stations. Figure 5(A) shows that the sensitivity of T_M to D_E is significantly and negatively correlated with incoming SW radiation at the end date of snow cover (R =− 0.29,P < 0.001). Secondly, albedo is the other factor that can feedback on temperature during snowmelt as discussed in previous studies (Groisman et al 1994, Bony et al 2006, Hall and Qu 2006, Fernandes et al 2009, Qu and Hall 2006, Flanner et al 2011). An earlier disappearance of snow influences temperature locally by diminishing surface albedo, causing more incoming SW radiation to be absorbed by the surface (Hall and Qu 2006). Here, we used albedo (black sky albedo (BSA)) from EOS-Terra-MODIS (MODerate-resolution Imaging Spectrometer) MCD43C3 version 5 products over the period 2001–2006 separated between snow free and snow cover (δBSA) to estimate the snow albedo change around the snowmelt (see section 2 and figure S1 available at stacks.iop.org/ERL/8/014008/mmedia). Figure 5(B) shows the relationship between the temporal sensitivity of T_M to D_E (see above) and δBSA across the 575 stations for which there is matching albedo data over the period 2001–2006 (R =− 0.44,P < 0.001). We found that the temporal sensitivity of T_M to D_E is significantly and negatively correlated with both incoming SW radiation and δBSA, verifying the larger difference in absorbed energy before and after snowmelt associated with larger snow albedo feedback on temperature (figures 5(A) and (B)). We took the product of incoming SW radiation by δBSA to present the extra absorbed energy by surface because of change in D_E (δASR), which exhibited a stronger relationship with the temporal sensitivity of T_M to D_E (R =− 0.50,P < 0.001) than that using incoming SW radiation or δBSA alone (R =− 0.29,P < 0.001;R =− 0.44,P < 0.001).

Zoom In Zoom Out Reset image size

In addition, the temporal sensitivity of T_O to D_O is much less correlated with incoming SW radiation during snow onset and δASR (R =− 0.21,P < 0.001;R = 0.06,P = 0.121; figure S2 available at stacks.iop.org/ERL/8/014008/mmedia). This may be related to the smaller incoming SW radiation during snow onset season (incoming SW radiation in 97% of stations smaller than 10 MJ d⁻¹) and δASR is limited (δASR in 97% of stations smaller than 2 MJ d⁻¹). However, the incoming SW radiation during snowmelt season (spring) exceeds 10 MJ d⁻¹ over most of the stations and δASR exceeds 2 MJ d⁻¹ over 51% of the stations. The albedo also has influence on the sensitivity of T_O to D_O, but the explanation is limited (R = 0.30,P < 0.001; figure S2 available at stacks.iop.org/ERL/8/014008/mmedia).

The average temporal sensitivity in temperate regions (south of 50°N) (−0.059 ± 0.055 °C d⁻¹, n = 358) is larger than that in boreal regions (north of 50°N) (−0.106 ± 0.075 °C d⁻¹, n = 217) (P < 0.001). This may be related to the fact that a higher incoming SW radiation could be available during snowmelt in boreal regions (16.1 ± 4.1 MJ d⁻¹) than in temperate regions (14.5 ± 3.8 MJ d⁻¹) over NH (P < 0.001), considering the later snowmelt date in higher latitude (figure S3 available at stacks.iop.org/ERL/8/014008/mmedia) (Choi et al 2010). In addition, δBSA in boreal regions (0.17 ± 0.10) is also higher than that in temperate regions (0.14 ± 0.1). This could explain why the magnitude of snow albedo feedback is larger in the high latitude region over NH (Serreze and Francis 2006, Chapin et al 2005, Euskirchen et al 2007, Bony et al 2006).

At Barrow station (latitude is about 70°N), Stone et al (2002) found that temperature in June increased by 0.9 °C when snowmelt date advanced by 14 d (table 2 in their paper). In other words, the sensitivity of temperature to D_E is −0.064 °C d⁻¹ at this site, which is a little smaller than the average result (−0.092 ± 0.054 °C d⁻¹) at 12 stations located between 69°N and 71°N in this study. Based on the seasonal cycles of satellite-derived albedo and air temperature, the sensitivity of surface albedo to temperature was estimated as −0.93 ± 0.06% K⁻¹ and −1.06 ± 0.08% K⁻¹ in Fernandes et al (2009) and Hall and Qu (2006) over the regions north of 30°N, respectively. In this study, the sensitivity of surface albedo to air temperature is estimated as the product of the sensitivity of D_E to temperature by δBSA of −0.90 ± 0.86% K⁻¹. In terms of radiative forcing, we used the product of sensitivity of D_E (D_O) to temperature by δASR as the magnitude of snow albedo feedback, which are 0.13 ± 0.13 W m⁻² K⁻¹ and 0.06 ± 0.09 W m⁻² K⁻¹ in spring and autumn, respectively. The total annual snow albedo feedback (spring plus autumn) over the land of Northern Hemisphere found in this study is 0.19 ± 0.17 W m⁻² K⁻¹, while the total cryosphere albedo feedbacks including snow and ice over NH was estimated to lie in the range of 0.3–1.1 W m⁻² K⁻¹ from satellite observation, and to average at 0.25 ± 0.17 W m⁻² K⁻¹ from 18 CMIP3 climate models (IPCC 2007, Flanner et al 2011). The estimate of our study, based on in situ records is comparable to these earlier estimates, although our results only include the local albedo feedback and are limited to the amount of in situ stations and their spatial distribution. The uncertainty in our estimate is larger, however, because it includes the variation across the 575 stations.

In this study, we used the extra absorbed energy at surface due to earlier snow disappearance (later snow onset) to explain the magnitude of the sensitivity of temperature to the end (onset) date of snow cover, assuming that the sensible heat flux allocation is related to absorbed energy at surface. Actually, earlier snowmelt can impact air temperature through changing the hydrological cycle. Snowmelt increases surface runoff in the beginning of the snowmelt period, but earlier snowmelt potentially leads to decreased soil water availability by the end of the snowmelt period (Hamlet et al 2005), particularly if frozen soils prevent snowmelt from recharging the soil moisture reservoir (Dunne runoff) (Dunne and Black 1970). As a consequence, an increase in the exchange of sensible heat between the earth surface and the atmosphere is expected; this will further contribute to an increase in T_M (Seneviratne et al 2006, Tchebakova et al 2002, Chapin et al 2008). In addition, we did not consider the cloud variation which could impact the solar radiation at the surface, and then impact the snow albedo feedback (Bony et al 2006, Qu and Hall 2006, Flanner et al 2011). More complex drivers such as sensible and latent heat flux allocation and cloud impacts on the snow albedo feedback should be included in future modeling studies.