Introduction

Snow cover over the Northern Hemisphere (NH) plays a crucial role in the Earth’s hydrology and surface energy balance and modulates feedbacks that control variations of global climate1. Previous studies found that snow cover in the NH has decreased significantly, i.e., the June snow cover extent (SCE) has decreased nearly twice as rapidly as the widely acknowledged September loss of sea ice extent from 1979 to 20112. This trend will continue in the future according to climate projections3,4. It coincides with north hemispheric warming and is indicative of a positive feedback of surface reflectivity on the climate3,5,6,7, which may affect the global biological and ecological systems further8,9.

In a warmer climate, the date of snow melt will advance in time10. Earlier snow melting can lead to major alterations in timing and volume of spring snowmelt runoff, with a possible increase in the incidence of catastrophic events such as spring flooding and summer droughts. Moreover, snow cover is tightly correlated with certain biological and ecological phenomena that depends on accumulated temperature, such as plant growth and animal migrations. During winter, snow cover is beneficial to agriculture by conserving the heat of the surface and thus, protecting the crops (mainly winter wheat) from the cold air. In high latitudes, warmer winter events damage the vegetation and reduce the following season crop yield11. Earlier snowmelt may reduce the availability of grass in its most nutritious form if animal migrations are not timed accordingly. In addition, changes in snow cover phenology may influence the seasonality in the terrestrial system through changes in the soil thaw and freeze dates12; as well as, in a warmer climate, result in permafrost degradation13.

Previous efforts have proved that snow cover phenology has remarkably changed in response to climate change. Reported effects at a local and regional scale include a shortening of snow duration days (Dd)14,15, an earlier snowmelt onset16 and an earlier snow end date (De)14,17. However, further work is still needed since published estimates of snow cover phenology mainly focused on the Northern pan-Arctic regions and conclusions were made from in-situ observations, model simulations or using a single source of images. Further studies are currently needed at a continental scale in order to better understand the changes in snow cover phenology, the mechanisms driving them and their impacts on the Earth’s climate system.

Recent studies have revealed significant changes in the land surface temperature induced by the vanishing cryosphere18,19,20,21, as well as the increased winter precipitation induced by changes in atmospheric circulation22 and human activity23,24. In particular, the near-surface at the NH high latitudes is warming at rates double of those at lower latitudes, due to the combined rapid loss of sea ice and snow cover in spring and summer20. In spite of the climate warming on average, an ostensibly large number of high-impact cold extremes have occurred in the mid-latitudes of NH over the past decade25. Snow cover phenology is highly sensitive to changes in temperature and precipitation. Therefore, understanding how these climate change events featuring the high spatial heterogeneity influence the snow cover phenology has meaningful consequences for water management, sustainable development of ecosystems and prediction of catastrophic climate related events. Thus, the objective of this study is to quantify and understand the spatial and temporal changes of snow cover phenology and to point out their causes and consequences. This knowledge is critical for assessing and projecting the future climate.

Five snow datasets are analyzed (Supplementary Table S1) in this study, including a reanalyzed dataset of daily snow depth generated by the Canada Meteorological Center (CMC)26, a binary daily snow cover mask derived from the Interactive Multi-sensor Snow and Ice Mapping System (IMS)27 and the Northern Hemisphere Weekly Snow Cover and Sea Ice Extent (NHSCE)27,28, 8-Day Level 3 snow cover fraction products (MOD10C2) derived from the Moderate Resolution Imaging Spectroradiometer Satellite (MODIS)29 and the snow water equivalent (SWE) derived from the Near-real-time Ice and Snow Extent (NISE) dataset30. A multi-data approach5 is employed in order to develop a combined snow cover phenology matrix that integrates snow cover phenology information detected from multiple sources of snow observations (see Methods, Supplementary Tables S2 and S3). Based on the seasonal cycle of snow cover over the NH31, we defined the snow cover accumulation season from the previous year (t−1) November to the current year (t) February and the snow cover melting season from March to June of current year (t). Regarding daily snow observations, the snow onset date (Do) is defined as the first five consecutive days on which snow was observed to cover the ground surface in accumulation season; De is defined as the last five consecutive days when snow cover was observed in the melting season. Dd is defined as the number of days from the onset date of snow cover to the end date of snow cover. The specific definitions of Do, De and Dd for each dataset are described in Methods.

Results

Observed changes in the snow cover phenology

Based on the NHSCE dataset, the NH presents no noticeable changes in Do, an earlier De at 95% confidence level (CL), as well as an overall shorter Dd between the winters of 1972/73 and 2007/0814. However, monthly SCE anomalies from 1972 to 2014 over the NH indicate that there is a significant reduction of SCE in summer and a notable increase of SCE in winter, especially after the year 2000 (Supplementary Fig. S1). Moreover, our results pointed out a contrast in the snow cover phenologies of middle and high latitudes over the NH in recent years. The snow onset date Do (Fig. 1a,d) over most parts of the NH high latitudes (40–70°N) was delayed by approximately 2.19 (±1.63) days except in latitudes lower than 40°N where Do advanced by 2.70 (±1.97) days. The most notable delay of Do (17.36 ± 6.67 days) occurred over Eastern Europe and Western Asia (the boxed area, Fig. 1a). Unlike the Do, the snow end date De (Fig. 1b,e) has advanced by about 9.66 (±2.35) days over high latitudes in Eurasia, Canada and the high elevation Tibet Plateau (TP) regions, but it was delayed by about 10.67 (±2.35) days in the mid-latitudes of Eastern Asia and Central North America (Fig. 1b). Combined spatial and temporal changes in Do and De led to an extension of snow cover duration days Dd (Fig. 1c,f) by over 10 days in latitudes below 40°N. At regions above 40°N, Dd shortened at a rate of 0.7 days decade−1, from 2001 to 2014.

Figure 1
figure 1

Changes of snow cover phenology over the NH from 2001 to 2014 derived from multi-data sets.

Changes of snow onset date Do (a,d), snow end date De (b,e) and snow duration days Dd (c,f) over the NH snow covered landmass, over 14 years. Changes are derived from the linear slope multiplied by the time span. Black dots in (a–c) indicate that the changes are significant at 90% CL. The zonal distribution in (d–f) are mapped at a 0.5 degree resolution in latitude. The error bars in (d–f) are calculated using equation (1) in Methods. The figure was created using ArcGIS42.

Full size image

De and Dd are strongly correlated (≥99% CL, Pearson correlation analysis) (r = 0.89) and the correlation between Do and Dd (r = 0.64) was also significant at a 95% CL. Attending to the small changes observed in Do (Supplementary Fig. S7), changes in Dd in the recent decades were mainly attributed to changes inDe.

In addition, yearly-averaged De decreased significantly (95% CL) over the NH snow covered landmass from 2001 to 2014 (Supplementary Fig. S7). Moreover, there were large differences in zonally averaged De anomalies between regions located in middle and high latitudes (Fig. 1b,e). The sign of the zonally averaged changes in De turns at around 52°N. The average change in De below 52°N was of 3.28 (±2.59) days, with maximum delays distributed around 35.5°N (9.78 ± 3.49 days). On the other hand, changes in De averaged at −5.11 (±2.20) days in the 52–75°N range, with maximum negative changes occurring at 72 °N (−8.59 ± 2.78 days). From 35.5°N poleward, De advances linearly with latitude.

Furthermore, the contrasting changes in De (Fig. 1b) observed between the NH middle and high latitudes are highly similar to the spatiotemporal distribution of changes in land surface albedo (as) observed during the melting season over the corresponding period ranging from 2001 to 2014 (Fig. 2a). The coefficients of linear correlation between the melting season De and as are above 0.8 over 91% of the studied area (Fig. 2b). In response to the earlier De observed in high latitudes during the melting season, as decreased due to lower land surface reflectivity, which resulted in the land-atmosphere system absorbing additional energy. On the other hand, delayed De observed in mid-latitudes during the melting season, resulted in opposite changes in as, which enhanced the cooling effects of the snow cover.

Figure 2
figure 2

Changes in land surface albedo during the melting season.

(a) Changes in surface albedo (%) during the melting season over 14 years and (b) correlation with De over the NH from 2001 to 2014. Black dots in (a) indicate that changes are significant at a 90% CL. The figure was created using ArcGIS42.

Full size image

Discussion

Attribution of snow cover phenology changes

The attribution analysis was carried out for the snow seasons from 2001 to 2013 limited by the availability of the Climatic Research Unit (CRU) Time Series (TS) dataset32. The analysis of linear correlations of snow onset date Do with land surface temperature (Ta) and precipitation (Pa) over the snow accumulation season (November-December-January-February) demonstrated that the Do was largely determined by Ta anomalies, pointing to a sensitivity of 0.28 days °C−1, as Pa was not statistically correlated with Do at a 95% CL (Supplementary Figs S8 and S9). The analysis of spatial patterns showed a general positive sensitivity of Do to Ta and negative sensitivity to Pa (Supplementary Fig. S10). Do was most significantly sensitive to Ta at around 40°N, as a result of the significant decrease in Ta in Eastern Europe and Western Asia (Fig. 3a). This spatial distribution is basically consistent with the observed long-term tendency of large-scale cooling trends of land surface temperature in winter over the mid-latitudes20,33. The decline of winter temperatures in mid-latitudes is attributed to changes in the Arctic system through changes in storm tracks, the Jet Stream, planetary waves and their associated energy propagation20. In addition, the recent Pacific Ocean cooling effect also contributed to a decrease of temperature in the mid-latitudes, since the tropical cooling effect on the extra-tropics is most pronounced in winter33.

Figure 3
figure 3

Attributions of snow cover phenology changes over the NH from 2001 to 2013.

Changes in (a) temperature during the accumulation temperature, Ta, (b) temperature during the melting season, Tm, (c) temperature over the entire snow season, Ts and (d) precipitation during the accumulation season, Pa, are derived from a linear regression model over the NH snow covered landmass from 2001 to 2013. (e) Yearly averaged snow end date De anomalies and contributions from Ta, Tm and Pa. (h) Zonally averaged De anomalies and contributions from Ta, Tm and Pa. Black dots in (a–d) indicate that changes are significant at a 90% CL. The figure was created using ArcGIS42.

Full size image

Changes in snow end date De were highly correlated with Ta, land surface temperature during melting season (March-April-May-June) (Tm) and Pa at a 95% CL, but Tm dominated (Fig. 3e). The De varies with Ta, Tm and Pa at −0.03 days °C−1, −0.31 days °C−1 and 0.10 days cm−1, respectively (Supplementary Fig. S11). The observed delays of De around 45°N were largely controlled by decreased Tm (Fig. 3b) as well as increased Pa (Fig. 3d) in this region, in which the magnitude of Pa was the main controlling factor of De (Fig. 3f). This pattern coincides with previous findings concerning large scale cold snaps, heavy snowfall and glacier events across the United States, Europe and East Asia25,34,35. Increased precipitation in the mid-latitude was attributed to afforestation24 and strengthening westerlies22,35. In addition, the advanced De around 70°N was largely explained by the poleward increase of Tm (Fig. 3f).

Snow duration days Dd anomalies were highly correlated in space and time with temperature over the entire snow season (Ts) and Pa at the 95% CL (Supplementary Figs S9 and S12). Both Ta and Tm significantly increase poleward above 60°N and decrease at different rates in the mid-latitudes (Fig. 3b, Supplementary Fig. S8b), which results in the observed contrasting Ts anomalies in middle and high latitudes (Fig. 3c). Contrasting Ts over the NH was mainly driven by Arctic amplification effects, as the warming magnitude at high latitudes was about two times as that of lower latitudes20,21. The rapid Arctic warming has contributed to dramatic melting of Arctic sea ice and spring snow cover, which is coincided with a period of ostensibly more frequent extreme weather events across the Northern Hemisphere mid-latitudes, including severe winters20. Dd varies with Ts and Pa at −0.85 days °C−1 and 0.14 days cm−1, respectively. The analysis of spatial patterns pointed out the overall negative sensitivity of Dd to Ts and its positive sensitivity to Pa (Supplementary Fig. S12). Even though changes in Dd are highly correlated with changes in Ts, correlations of Dd with Tm and Ta are not statistically significant.

Radiative forcing and albedo feedbacks of snow cover phenology changes

Land surface change is expected to contribute substantially to warming trends in the NH high latitudes36. To quantify the snow radiative forcing (SnRF) and albedo feedback induced by snow cover phenology changes, both the radiative kernel method37,38,39 and the Clouds and Earth’s Radiant Energy System (CERES) observations at the top of atmosphere (TOA) were employed in this study. The radiative kernel approach allows to separate the radiative response to different climate parameters and to decompose the feedback into radiative and climatic response components. In this study, the albedo radiative kernel was used to evaluate the instantaneous perturbation to the Earth’s shortwave (SW) anomaly induced by snow cover phenology changes. The CERES observations were used to study the snow cover phenology anomaly induced radiative forcing in total TOA forcing anomalies.

Applying the radiative kernel approach with equation (2), we inferred that changes in SnRF induced by the snow season variability from 2001 to 2013 were recorded to be 0.16 (±0.004) Wm−2 over the NH. Distinct anomalies are observed at middle and high latitudes, where changes in the snow radiative forcing during the accumulation season (SnRFa) was of 0.01 (±0.001) Wm−2 and changes during melting season (SnRFm) was of 0.31 (±0.011) Wm−2. Compared with changes in SnRFa, those in SnRFm are the main contributor to changes in the SnRF in both space and time, (Fig. 4d,e). As displayed in Fig. 4c, notable positive anomalies of SnRF were observed on the landmass near the Arctic Ocean, at altitude on the Rocky Mountains and in the TP regions and negative changes were observed in the mid-latitudes of Eastern Asia and Central United States which is highly identical with spatiotemporal distribution of SnRFm displayed in Fig. 4b.

Figure 4
figure 4

Spatiotemporal distribution of snow radiative forcing SnRF (in Wm−2) and contributions from the accumulation and the melting seasons from 2001 to 2013 using the radiative kernel approach.

Estimates of SnRF induced by (a) changes in the snow onset date Do, (b) the snow end date De and (c) the snow duration days Dd from 2001 to 2013. (d) SnRF and contributions from the accumulation season SnRFa and the melting season SnRFm. (e) Zonally averaged SnRF and contributions from the SnRFa and SnRFm. Value in (a–c) were averaged from two SnRFa, SnRFm and SnRF results obtained using the Community Atmosphere Model (CAM3) and the Atmosphere Model 2 (AM2), respectively. Black dots in (a–c) indicate that changes are significant at a 90% CL. The figure was created using ArcGIS42.

Full size image

Direct estimations using CERES observations further confirmed that changes in De largely contributed to TOA SW anomalies (Fig. 5a and Supplementary Fig. S13) as suggested by the correlation coefficient of 0.87 between De and the TOA SW flux anomaly during the melting season, from 2001 to 2013. Changes in De influenced the SW flux more than the longwave (LW) flux (Fig. 5b). This resulted in the positive anomaly of TOA net flux and led to further warming. The De-induced SnRFm of 0.31 (±0.01) Wm−2 accounted for 51% of the total TOA SW flux anomaly (0.61 Wm−2) and 63% of the TOA Net flux changes (0.50 Wm−2) in melting season over the NH snow-covered landmass from 2001 to 2013. In addition, the observed TOA SW flux changes are highly consistent with the observed De anomalies from mid-to-high latitudes, as a weakened cooling effect is observed in the high latitudes as well as an enhanced cooling effect in the mid-latitudes, which indicates that changes in De may further influence and control local and regional climate change.

Figure 5
figure 5

Linear correlation of the snow end date De anomalies with the TOA (a) SW flux, (b) LW flux and (c) Net flux over the NH snow covered landmass from 2001 to 2013.

SW, LW and Net flux were derived from CERES observations. Linear correlations in (a–c) are significant at the 95% CL.

Full size image

NH landmass warming from 2001 to 2013 was estimated from the Goddard Institute for Space Studies (GISS) Surface Temperature Analysis40 at 0.07 °C. Combining the amplitudes of SnRF, SnRFa and SnRFm with this warming estimate (ΔT) yields a NH snow albedo feedback (SnRF/ΔT) of 2.29 (±0.06), 0.14 (±0.01) and 4.43 (±0.14) Wm−2K−1, respectively. The NH Cryosphere radiative forcing caused by snow cover variability between 1979 and 20086 was estimated at 0.27 (0.11–0.48) Wm−2 with a warming of 0.79 °C from GISS Surface Temperature Analysis, which yields a snow albedo feedback of 0.34 (0.14–0.61) Wm−2K−1. Compared with results from Flanner et al.6, our results indicate that even though the NH snow radiative forcing SnRF from 2001 to 2013 (0.16 Wm−2) was about 60% of SnRF from 1979 to 2008 (0.27 Wm−2 Wm−2), the NH snow albedo feedback during 2001–2013 (2.29 Wm−2K−1) was about six times of those from 1979 to 2008 (0.34 Wm−2K−1) due to a lower warming extent.

Our results based on the snow cover phenology matrix obtained from several datasets from 2001 to 2014 indicate contrasting snow cover phenology characteristics in the NH middle and high latitudes. These results can help to better understand the spatial and temporal changes of snow cover phenology at a continental scale, in the context of the recent climate change and help to predict it. Climate models made at a continentental scale should take the contrasting features of snow cover phenology into consideration to reduce uncertainties in climate projections. However, the availability of long, high-resolution time-series of snow cover observations from satellites is expected to improve the future study of snow cover phenology. The warming trend is very likely to continue with the accumulation of greenhouse gas in the atmosphere; as well as the intensified effects of the Arctic Amplification resulting from the vanishing cryosphere snow and sea ice. Accordingly, the contrasting changes in the snow cover phenology are very likely to continue in the near future. Moreover, snow cover phenology is largely determined by temperature and precipitation anomalies. According to the similarities between the distribution of snow radiative forcing, the observed TOA SW flux anomalies and the observed contrasting anomalies of De over the NH, it is necessary to further investigate how these energy budget anomalies control climate change.

Methods

Study area

In order to discern variations in the NH snow cover phenology, the snow cover regions (grid cells) was defined as the area ranging between 32°N and 75°N, excluding grid cells of permanent snow cover. The northernmost part of Greenland was excluded from the analysis as there are few products providing reliable snow cover information due to the complex coastal topography and the difficulty in discriminating snow from ice5. The region below 32°N was also excluded from the analysis as patchy snow or thin snow covering vegetated surfaces may be missed in the Normalized Difference Snow Index (NDSI) calculated using an algorithm for snow mapping that processes visible and near-infrared images, such as MODIS41. In addition, snow can melt and reappear over the course of the snow season. At higher latitudes or altitude, this may only occur at the beginning or at the end of the season, whereas at lower latitudes this may occur throughout the winter14. To allow for inter-annual comparisons of snow cover phenology, this study is restricted to stable snow-covered regions where snow covered the ground for at least 60 days each year from 2001 to 2014.

Snow cover phenology retrieval

Five snow datasets derived from satellite observation and reanalyzes are analyzed (Table S1) in this study. For daily CMC snow depth observations, Do is defined as the first five consecutive days and De is defined as the last five consecutive days when snow depth exceeds 1 cm, over a given year (t). For the daily IMS binary SCE mask dataset, Do and De are defined as the first and last five consecutive images when pixels were marked as 1 in the records, over a given year (t). For the daily NISE SWE dataset, Do and De are defined as the first and last five consecutive days when the snow depth exceeds 2.5 cm, over a given year (t). For the weekly NHSCE binary SCE mask dataset, we first identified the date range (i to i + 6) of the first frame (n) when snow cover occurred and disappeared, respectively. Then, over a given year (t), Do and De are defined at i + 3, relative to the first frame when snow cover occurred and disappeared, respectively. For the 8-day MOD10C2 snow cover fraction dataset, we first identified the interval (i to i + 7) limits between the frames when the snow cover fraction is superior to 0 and subsequently equals 0, then defined Do and De at i + 3.5 days of first frame when the snow cover fraction is superior to 0 and equal to 0, respectively.

In order to keep all the information held in the five individual datasets as well as to establish spatial comparisons and causal analyses between snow phenology information obtained from each one of them, we first identify Do, De and Dd from individual datasets without regridding. Then, the results obtained from individual datasets are regridded at 0.5 degree spatial resolution using a resampling method of “average” with the help of gdalwarp (http://www.gdal.org/gdalwarp.html). The applied resampling method computes the average of all non-NODATA contributing pixels in the domain of our study.

Multi-data snow cover phenology retrieval and uncertainty analysis

A multi-data approach5 is employed to develop an integrated snow cover phenology matrix and reduce the uncertainty from individual datasets. Taking De as an example, individual datasets of snow observations suggest different values and spatiotemporal distributions of De (Supplementary Figs S1 and S2, supplementary Tables S2 and S3) depending on the spatial resolution, the method (or algorithm) used to detect snow cover and the different definitions of snow cover. In order to achieve internal consistency, De estimates derived from each individual dataset were converted to standardized z-score anomalies using the mean and standard deviation values observed over a given period (Supplementary Fig. S4). The consistency of each dataset was evaluated by computing the correlation and the root-mean-square error (RMSE) of the multi-dataset mean, excluding the dataset being verified (Supplementary Tables S4 and S5). The final integrated multi-dataset De series was obtained by averaging the anomaly series from 2001 to 2014. The final De series was then converted back to a De value for each day of the year using the mean and standard deviation obtained from the MODIS dataset since MOD10C2 constitutes a source of consistent and objective snow cover phenology estimates, derived from high resolution optical satellite data, compared with the CMC, NHSCE, IMS and NISE snow datasets used in this study. This approach was also applied at pixel level to help remove the individual De series with poor data quality from the final averaged anomaly series. The correlation and RMSE of each De anomaly series with the average anomaly series from the other four snow data sets are displayed in supplementary Fig. S5. An estimate of the uncertainty in De is obtained from the standard error (SE),

which depends on the standard deviation s of the n data sets included in the average anomaly. The multi-data Do and Dd series were estimated following the same process.

For interannual and spatiotemporal anomalies in Do, De and Dd, we first calculated interannual and spatiotemporal changes from each individual dataset. Then, averaged the individual changes with significance above 90% CL to obtain the final multi-data anomalies. The estimate of the uncertainty in the final multi-data changes is obtained from SE (equation 1), which depends on the standard deviation s of the n data sets included in the final multi-data changes calculation. In addition, we apply 2 s control limits (−2 and +2 times standard error) in the calculation of multi-data anomalies. Individual anomalies beyond this range were also excluded in this study.

Snow radiative forcing (SnRF) calculation using albedo radiative kernels

The albedo radiative kernel is expressed as the TOA net shortwave anomalies associated with a 1% change in the land surface albedo as. In our study, the albedo radiative kernel — estimated using radiative transfer algorithms from the Community Atmosphere Model (CAM3) of the National Center for Atmosphere Research as well as the Atmosphere Model 2 (AM2) of the Geographical Fluid Dynamic Laboratory, developed by Shell et al.38 and Soden et al.37— is used to quantify the SnRF. We define the SnRF as the instantaneous perturbation to Earth’s TOA SW anomalies consequent of the snow season variability, which can be quantified as as anomalies driven by the disappearance of snow cover during the melting season and by the presence of snow cover during the accumulation season. Thus, the time (t) dependent SnRF within the study area of area A, which is composed of gridcells r, can be represented as follows:

Here, S is the snow cover fraction over the study area. ∂αs/∂S is the rate of variation in surface albedo with the snow cover changes. ∂F/∂αs is the response of TOA SW flux variation to as changes. We assume that (monthly and spatially varying) ∂αs/∂S and ∂F/∂αs are constant with snow cover fraction S and surface albedo αs, respectively. Then, ∂αs/∂S can be replaced with mean albedo contrast induced by snow cover anomaly and ∂F/∂αs can be obtained from the albedo radiative kernels6,39.

Additional Information

How to cite this article: Chen, X. et al. Observed contrast changes in snow cover phenology in northern middle and high latitudes from 2001–2014. Sci. Rep. 5, 16820; doi: 10.1038/srep16820 (2015).