Satellite observed changes in the Northern Hemisphere snow cover phenology and the associated radiative forcing and feedback between 1982 and 2013

Quantifying continental-scale changes in snow cover phenology (SCP) and evaluating their associated radiative forcing and feedback is essential for meteorological, hydrological, ecological, and societal purposes. However, the current SCP research is inadequate because few published studies have explored the long-term changes in SCP, as well as their associated radiative forcing and feedback in the context of global warming. Based on satellite-observed snow cover extent (SCE) and land surface albedo datasets, and using a radiative kernel modeling method, this study quantified changes in SCP and the associated radiative forcing and feedback over the Northern Hemisphere (NH) snow-covered landmass from 1982 to 2013. The monthly SCE anomaly over the NH displayed a significant decreasing trend from May to August (−0.89 × 106 km2 decade−1), while an increasing trend from November to February (0.65 × 106 km2 decade−1) over that period. The changes in SCE resulted in corresponding anomalies in SCP. The snow onset date (Do) moved forward slightly, but the snow end date (De) advanced significantly at the rate of 1.91 days decade−1, with a 73% contribution from decreased SCE in Eurasia (EU). The anomalies in De resulted in a weakened snow radiative forcing of 0.12 (±0.003) W m−2 and feedback of 0.21 (±0.005) W m−2 K−1, in melting season, over the NH, from 1982 to 2013. Compared with the SCP changes in EU, the SCP anomalies in North America were relatively stable because of the clearly contrasting De anomalies between the mid- and high latitudes in this region.


Introduction
The snow cover extent (SCE) and snow cover phenology (SCP) over the Northern Hemisphere (NH) play a crucial role in the Earth's hydrology and surface energy balance, and modulate the feedbacks that control variations in global climate (Frei et al 2012, Peng et al 2013, Brands 2014, Chen et al 2015. The Eurasian snow cover in October is closely related with the December-February mean boreal winter climate via the 10 m wind speed and significant wave heights in the North Atlantic and adjacent seas (Brands 2014).
Changes in the Eurasian fall snow cover could be a potential driver of changes in the Arctic Oscillation. Such changes may also provide a link between the decline in sea-ice in the Arctic during summer and the unusual atmospheric circulation the following winter (Wegmann et al 2015). Trends in the duration or extent of the snow cover are expected to exert feedback on the temperature trends (Peng et al 2013). In a warmer climate, snow will melt earlier in the year; and this has already been occurring in some places (Barnett et al 2005). Moreover, snow cover changes are also likely to impact the carbon pool stored in the permafrost (Gouttevin et al 2012). SCP is highly correlated with various biological phenomena that depend on the accumulated temperature, such as animal migrations and vegetation growth (Bokhorst et al 2009). In addition, the surface radiative and turbulent fluxes in spring and summer have strong correlations with satellite-observed SCE (Barnett et al 2005). Changes in the snow cover in the NH landmasses enhance climatic anomalies through radiative forcing and snow-albedo feedback (SAF) effects. As surface air temperature increases, the snow cover over the NH extra-tropical region retreats, which causes the land surface to reflect less solar radiation. The extra radiation is absorbed, resulting in enhanced warming, especially in the regions undergoing snow and ice reduction (Hall 2004). This SAF accounts for approximately half of the additional net incoming solar radiation associated with the retreat of the NH cryosphere in equilibrium climate change simulations (Hall 2004, Qu and. Thus, monitoring the SCE and SCP is essential to assess the impacts of climate change and projections of future climate, especially in the snow-dominated regions. Previous studies have proven that SCE in the NH has experienced a well-documented rapid decrease since satellite observations began in the late 1960 s. Climate projections suggest that the SCE will continue shrinking in the future Robinson 2011, IPCC 2013), coinciding with the hemispheric warming and indicating the positive feedback of surface reflectivity on climate (e.g. Brown et al 2010, Brown and Robinson 2011, Flanner et al 2011. In response to the variability of SCE, the SCP has also experienced remarkable changes. Such changes include the shortening of snow cover duration days (D d ) (e.g. Whetton et al 1996, Choi et al 2010, with an earlier snow melt onset , an advanced snow end date (D e ), and a delayed snow onset date (D o ) (Choi et al 2010), at both the local and regional scales. However, only a few of these studies have explored the effects of long-term SCP changes on the radiative forcing and its associated feedback. The effects of longterm SCP changes on the Earth's climate system and energy budget are still unclear. Additionally, most published studies are based on in situ observations or models in which the in situ measurements are highly dependent on a particular location (latitude and elevation), hence limiting the model outputs with assumptions. Moreover, published estimates of SCP focus mainly on high-latitude regions of the NH, and conditions at the hemispheric scale are rarely documented.
Given the importance of SCP in the Earth's climate system and in the current climate change background, it is necessary to quantify and understand its spatial and temporal anomalies, and to clarify its associated radiative forcing and feedback onto the Earth's climate system, using long-term satellite observations. Therefore, the objective of this study is to quantify long-term changes in SCP based on satellite observations and to explore its associated snow radiative forcing (S n RF) and albedo feedback over the NH snow-covered landmass between 1982 and 2013. In order to achieve this objective, we first quantified the SCE anomalies over the NH. Then, we retrieved the D o , D e , and D d based on the observed changes in SCE. Finally, we calculated the associated S n RF and albedo feedback using the radiative kernel technology.

Datasets
To investigate the latest SCE changes over the NH, we employed the widely used monthly SCE chart data generated by the Rutgers University Global Snow Lab and the binary snow cover mask data derived from the Climate Data Record of Northern Hemisphere Snow Cover Extent (NHSCE) at a spatial resolution of 25 km in this study. Both the chart data and the binary snow cover mask data were developed by Robinson et al (1993).
The NHSCE data was derived from manual interpretation of Advanced Very High Resolution Radiometer (AVHRR), Geostationary Operational Environmental Satellite, and other visible-band satellite data (Helfrich et al 2007). An uncertainty analysis for NHSCE in March and April was conducted, estimated from an inter-comparison of multiple datasets by Brown and Robinson (2011). It indicated a 95% confidence level in SCE estimations of ±3-5% over the satellite era begun in the late 1960 s. Since the binary value in the NHSCE dataset indicates a 50% or larger probability of occurrence of snow in each pixel, this data is more appropriate for large-scale snow cover studies. However, this dataset is still valuable in large-scale snow related studies. For instance, Choi et al (2010) derived the information on snow seasons from NHSCE and found that the duration of the average NH full snow season has decreased significantly between the winters of 1972/73 and 2007/08; Derksen and Brown (2012)  To calculate the albedo contrast (Δa s ) induced by changes in SCP, we employed the newly released longterm, high-quality, and spatially complete Global LAnd Surface Satellite (GLASS) land surface albedo (a s ) product  in 0.05°spatial resolution over the 1982-2013 period, and the MODerate resolution Imaging Spectroradiometer (MODIS)derived global gap-filled and snow-free a s product (MCD43GF) (Moody et al 2008) in about 0.0083°spatial resolution between 2001 and 2013. The GLASS a s is a real time series of surface albedos and the MCD43GF a s is a snow-free albedo. We then obtained in the present study the Δa s time series (yearly and monthly) from the GLASS a s , minus the climatology snow-free a s obtained from MCD43GF. The GLASS a s dataset was produced from both AVHRR and MODIS data  and has previously been used to quantify the radiative forcing of snow melt over Greenland  and snow cover changes over China (Chen et al 2016), with high quality and fine spatial resolution. The spatially complete MCD43GF snow-free a s was produced using an ecosystem-reliant temporal interpolation technique that retrieves missing data with an error of 3%-8% (Moody et al 2008).
In order to quantify the S n RF and the feedback caused by a s changes induced by SCP anomalies, the albedo radiative kernel expressed as the net shortwave (SW) anomalies at the top of the atmosphere (TOA) and associated with a 1% change in the a s was used. Such data was estimated using radiative transfer codes from the Community Atmosphere Model (CAM3) of the National Center for Atmosphere Research and the Atmosphere Model 2 (AM2) of the Geographical Fluid Dynamic Laboratory developed by Shell et al (2008) and Soden et al (2008).

SCP retrieval
To discern variations in the NH SCP, the snow cover region (grid cells) was defined as the area between 30°N and 75°N, excluding grid cells of permanent snow cover. Greenland was excluded from the analysis. There are indeed few products providing reliable snow cover information, due to the complex coastal topography and the difficulty of discriminating snow from ice (Brown et al 2010). Latitudes below 30°N were also excluded from the analysis because patchy snow or thin snow cover on vegetated surfaces may be missed in visible and near-infrared images (Riggs et al 2006). In order to explore the detailed information about changes in SCP over the NH, we analyzed changes in Eurasia (EU) and North America (NA) separately. Considering the fluctuation of SCE over the years, we confined our study area to stable snow-covered regions where snow covered the ground for at least 8 weeks (about 60 days) in 75% of the years between 1982 and 2013, following Choi et al (2010). The study area is displayed in figure 1.
In order to compare the SCP between different years, we used the hydrological year in this study. Based on the seasonal cycles of observed snow cover for northern EU and NA (Brutel-Vuilmet et al 2013), we defined the hydrological year of the NH as the period between September of a given year (t−1) and August of the following in year (t). Moreover, we defined the accumulation season of snow cover as extending from September of a given year (t−1) until the next February of the following year (t), while the melting season ranged from March to August of a same year (t). For the binary and weekly NHSCE dataset, we first identified the date range (i-i+6) of the first frame (m) when the snow cover first appeared in the accumulation season, as well as the date range ( j-j+6) of the last frame (n) when the snow cover first disappeared in the melting season. Then, D o and D e of the year (t) were defined as the date i+3 of the frame (m) in the accumulation season and j+3 of the frame (n) in the melting season. D d was defined as the time span between D o and D e .

Contribution analysis
Changes in the SCP over the NH are determined using the SCP anomalies over the EU and NA. To compare the contributions of SCP anomalies in EU and NA consistently, SCP variations over the NH, EU, and NA were converted into the standardized SCP anomalies (z-score) in the contribution analysis using the mean and the standard deviation of each series between 1982 and 2013. The contributions of EU and NA SCP anomalies to the NH SCP variability were computed annually to display the changing influence of EU and NA SCP anomalies on the underlying NH SCP variability.
The contributions of EU and NA SCP anomalies to the NH SCP variability were computed by regressing the annual series of NH SCP z-scores against the z-scores time series of the EU and NA SCP series. The resulting regression coefficients were then multiplied by the time series of EU and NA SCP z-scores to derive the contributions of the EU and NA SCP to the NH SCP variability. This method was employed by Pederson et al (2013) to identify regional patterns and proximal causes of the recent snowpack decline in the Rocky Mountains.

Radiative forcing and feedback calculation
The radiative kernel method is widely used to quantify the radiative forcing induced by a s changes. For example, Flanner et al (2011) quantified the radiative forcing from albedo change induced by changes in the NH cryosphere. Chen et al (2015) calculated the SCP change and its associated radiative cooling. In the present study, the albedo radiative kernels developed by Shell et al (2008) and Soden et al (2008) were used to quantify the S n RF and the feedback caused by the albedo change, which is induced by the SCP anomaly. Compared with analytical models, the radiative kernel methods are more accurate because they capture more effectively the functional dependence of the planetary albedo upon the surface albedo (Qu and Hall 2014), which is one of the most important parameters in the quantification of the albedo-induced radiative forcing and feedback.
We defined the S n RF as the instantaneous perturbation to Earth's TOA net SW anomalies consequent to the snow season variability. The a s anomalies driven by the disappearance of snow cover during the melting season and by the presence of snow cover during the accumulation season allow to quantify this snow season variability. Thus, the time (t)-dependent S n RF (W m −2 ), within a region R (here, the NH snow-covered landmass) of an area A composed of grid cells r, can be represented as follows: R n s s where S is the snow cover fraction(SCF) over the study area, ∂α s /∂S is the rate of variation of the surface albedo with the snow cover change, and ∂F/∂α s is the response of the TOA net SW anomalies to a s change.
We assumed that the monthly and spatially varying ∂α s /∂S and ∂F/∂α s are constant for the SCF S and the surface albedo α s , respectively. Then, ∂α s /∂S can be replaced with the mean albedo contrast induced by the snow cover anomaly and ∂F/∂α s can be obtained from the albedo radiative kernels (Flanner et al 2011, Chen et al 2015. The strength of the albedo feedback induced by changes in SCP can be quantified using the S n RF/T s values. In our study, S n RF/T s is defined as the amount of additional TOA net SW radiation averaged over the decreases in NH extra-tropical landmass a s , in association with a 1°C increase in surface air temperature (T s ) (Hall andHall 2007). We first calculated the S n RF for each month, then obtained the annual-mean S n RF by averaging monthly S n RF values from September to August, as defined in section 2.2.1. Similarly, we obtained the mean S n RF per accumulation season by averaging S n RF values from September to February and the mean S n RF per melting season by averaging S n RF between March and August.

Observed changes in SCE
The monthly SCE anomalies over the NH between 1982 and 2013, including EU and NA, are shown in figure 2. Detailed changes of monthly SCE over the NH are presented in table 1. As shown in figure 2 and table 1, the May-June-July-August (MJJA) period experienced a significant SCE reduction (−0.89 10 6 km 2 decade −1 ) at the 95% significance level, especially after the year 2008 (−2.55 10 6 km 2 decade −1 ). However, the monthly SCE variability over the 1998-2008 period in MJJA presented substantially different temporal anomalies in EU and NA. EU experienced a continuous SCE reduction after 1998 (−0.75 10 6 km 2 decade −1 ), contrasting with the SCE increase in NA (0.03 10 6 km 2 decade −1 ).
In contrast with the marked MJJA SCE reduction, the SCE increased over November-December-January-February (0.65 10 6 km 2 decade −1 ), especially after the year 2002. This temporal pattern in the accumulation season is consistent with the pattern of  Seasonal SCF anomalies were used to further investigate the spatial pattern of SCE anomalies over the NH from 1982 to 2013. For each pixel, the SCF was expressed as the percentage of snow cover appearance in a given season. The temporal anomalies of SCE over the NH were further verified by the spatial anomalies of SCF displayed in figure 3. The spring SCF in EU and Western NA, as well as the summer SCF over the entire NH, significantly decreased over the 1982-2013 period, during which there was a reduction of SCF in summer reaching up to 50% in high latitudes around 70°N at the 95% significance level. In contrast, the SCF in autumn and winter experienced a significant increase with magnitudes ranging 5%-30% at all latitudes over the NH. Increased SCF, especially in EU, extends the findings of Bulygina et al (2011) to a larger spatial scale, showing that the Russian territory has been experiencing an increase in snow depth since 1966, in the context of global temperature rise and sea ice reduction in the northern hemisphere. Moreover, the spatiotemporal distribution of SCF in autumn and winter coincides with the increasing snow cover and widespread boreal winter cooling observed over the last two decades (Cohen et al 2012). The 32 year averaged D e displayed a mean DOY value of 132.66 over the NH stable snow-covered landmass. The values for EU and NA were 128.74 and 139.86, respectively (supplementary figure s1(c)). The 32 year changes in D e over NH, EU, and NA were −6.12, −7.17 and −3.78 days decade −1 (no statistical significance), respectively. As shown in figure 4(c), almost all of the pixels over the NH snow-covered landmass, from low to high latitude, had an earlier D e over the 1982-2013 period. However, the contribution analysis indicated that D e changes in EU accounted for 73% of the D e anomaly over the NH from 1982 to 2013, while D e changes in NA only accounted for 27% of the NH D e changes ( figure 5(b)). This result extends the findings of Peng et al (2013) from the level of in situ observation to that of large continuous spatial coverage. Peng et al (2013) found, unlike their early D e in EU (−2.6±5.6 days decade −1 ), that the trends in D e over the NA (0.1±5.8 days decade −1 ) exhibited a fragmented spatial pattern from 1978 to 2006.

Changes in SCP over the NH
The combined spatial and temporal changes in D o and D e led to longer D d in central NA, Western Europe, and Eastern Asia ( figure 4(e)). The 32 year averaged D d displayed a mean value of 196.53, 191.11, and 206.48 days over the NH, EU, and NA, with a shortening trend between 1982 and 2013 of 1.04, 1.17, and 1.33 days decade −1 , respectively. Although a statistical analysis based on the linear regression model revealed that there were no significant changes in D d from 1982 to 2013, the inter-annual variability of D d still presented a sharp shortening trend after the year 2006, which was caused by the delayed D e over this period. These findings largely coincide with the significant reductions of SCE in May and June SCE in Derksen and Brown (2012). The high contributions from EU in the D o and D e anomalies result in a correspondingly high D d anomaly in EU. The contribution analysis indicated that EU accounted for 76% of the D d changes over the NH, against only 24% for NA over this period ( figure 5(c)).
The correlation analysis of D o , D e , and D d reveal significant relationships (99%, Pearson correlation analysis) between D e and D d (r=0.78, 0.70, 0.79 for the NH, EU, and NA, respectively) over the 1982-2013 period. The correlations between D o and D d (r=0.69, 0.73, 0.63 for the NH, EU, and NA, respectively) are also significant at the 95% level. Considering that there are no statistically significant changes in D o , the SCP changes over the NH from 1982 to 2013 are mostly attributed to the changes in D e , as demonstrated by future snow cover changes in EU.

Changes in the snow radiative forcing S n RF and the albedo feedback
The spatial distributions of annually averaged clearsky S n RF, accumulation season S n RF a and melting season S n RF m and their changes over the NH snowcovered landmass between 1982 and 2013 are shown in supplementary figure s2 and figure 6. We calculated the S n RF per month, and averaged the monthly S n RF values from the accumulation season and from the  melting season to obtain the value of S n RF a and S n RF m , respectively. Determined by the spatial pattern of SCP anomalies and the magnitude of the albedo radiative kernels in the high-to-low latitudes, the S n RF is generally larger over the high latitudes in the melting season ( figure s2(b)), as well as in the entire snow season (figure s2(c). This is a consequence of the large Δa s induced by a larger SCF compared with the lowto-mid latitudes, which shows a different distribution from the S n RF in the accumulation season ( figure s2(a)).
The spatial distribution of S n RF m changes (figure 6(c)) presents weakened cooling effects with smaller magnitude of S n RF m in high latitude landmass near the Arctic Ocean, and strengthened cooling effects with larger magnitude of S n RF m in the central United States from 1982 to 2013. This coincides with the observed spatiotemporal anomalies of D e (figure 4(c)). However, the S n RF m changes in NA are not statistically significant, which may be caused by the contrast between S n RF m changes over the midand high latitudes in this region. The spatial distribution of S n RF change is displayed in figure 6(e), which shows that longer D d in the central United States resulted in a correspondingly strengthened S n RF in this region. However, the longer D d distributed over Europe and Eastern Asia does not match well the distribution of S n RF in this region.
The 32 year interannual variations of clear-sky S n RF a , S n RF m , and S n RF from 1982 to 2013 are shown in figures 6(b), (d) and (f). The 32 year changes in S n RF a , S n RF m , and S n RF over the NH are 0.09 (±0.009), 0.12 (±0.003), and 0.10 (±0.005) W m −2 , respectively, at the 95% significant level. The statistical analysis revealed that EU accounted for 68%, 74%, and 70% of the changes in S n RF a , S n RF m , and S n RF,  The ranges of S n RF a , S n RF m , and S n RF combined with this warming yield changes in the accumulation season, melting season and entire snow season albedo feedback (S n RF/T s ) of 0.18 (±0.008), 0.21 (±0.005), and 0.17 (±0.008) W m −2 K −1 , respectively, over the NH.

Discussion and conclusion
SCP is one of the most sensitive factors to changes in climate. In this study, D o , D e , and D d were retrieved, based on observed NHSCE snow cover data. The associated radiative forcing and feedback were quantified partly using the observed surface albedo dataset and the albedo radiative kernels technology.
Compared with the few changes in D o , the 32 year anomalies in D e over the NH displayed a marked trend of −1.91 days decade −1 at the 95% significance level. The contribution analysis indicated that the D e change in EU (−2.24 days decade −1 ) accounted for 73% of the D e changes over the NH. The combined spatial and temporal changes in D o and D e led to a shortening trend of 1.04 days decade −1 over the NH from 1982 to 2013. However, earlier D o still resulted in longer D d in central NA, Western Europe, and Eastern Asia.
Driven by the SCP anomalies over the NH, the associated radiative forcing in the accumulation season, melting season, and entire snow season from 1982 to 2013 are 0.09 (±0.009), 0.12 (±0.003), and 0.10 (±0.005) W m −2 , respectively. Even through there are no statistically significant changes in D o , the anomaly in accumulation season radiative forcing is still significant at the 95% level. The statistical analysis revealed that EU accounts for 68%, 74%, and 70% of the changes in S n RF a , S n RF m , and S n RF over the NH between 1982 and 2013, which is consistent with the high contributions from EU to the SCP changes. However, S n RF m changes in NA are not statistically significant, due to the contrasting S n RF m changes in the mid-and high latitudes induced by the distribution of D e anomalies in this region.
Compared with previous research, the present study investigated the SCP changes, quantified the associated radiative forcing and feedback, and explored the contributions from EU and NA, which is important for current climate change studies and future climate projections. However, several unresolved issues remain in our study. S n RF is influenced by several factors, including the a s change induced by snow, the local insolation, and the atmospheric state (Flanner et al 2011). Limited by the spatial resolution of NHSCE snow cover dataset, the snow range used in this study includes variability caused by both unresolved snow cover variability and variability in the influence of vegetation on albedo contrast.
Considering the rapid responses of the permafrost and vegetation to experimentally increased snow cover (Johansson et al 2013), a separation of vegetation effects from S n RF quantification is necessary in future studies. In addition, since the NH warming trend is very likely to continue due to increasing greenhouse gas concentrations and intensified Arctic amplification effects, the SCP is also likely to experience rapid changes in the near future, which should be a high concern in climate change projections.