Distribution, attribution, and radiative forcing of snow cover changes over China from 1982 to 2013

Snow is an integral component of the global climate system. Owing to its high albedo and thermal and water storage properties, snow generates important linkages and feedbacks through its influence on surface energy and moisture fluxes, clouds, precipitation, hydrology, and atmospheric circulation. By virtue of its radiative and thermal properties that modulate transfers of energy and mass at the surface–atmosphere interface, snow affects the overlying atmosphere and therefore plays an important role in the complex web of feedbacks that control local to global climate (Frei et al. 2012). These functions could be significantly affected if snow cover is reduced in the future, as is an expected result of global warming through enhancement of the greenhouse effect (Whetton et al. 1996).

The impact of global climate change is the main topic of climatological research worldwide (Jaak 1997). The Intergovernmental Panel on Climate Change (IPCC) has reported that the snow cover extent (SCE) has obviously declined in the Northern Hemisphere (NH), and it has predicted that anthropogenic global warming will be greatest in northern snow-covered areas (IPCC 2013). Both observations (Brown and Robinson 2011; Derksen and Brown 2012) and climate models (Derksen and Brown 2012; McCabe and Wolock 2009; Roesch and Roeckner 2006) suggest an apparent decline in spring SCE and a strong direct influence of snow cover on the surface–atmosphere energy budget, particularly in the radiative forcing at the surface and the top of atmosphere (TOA). When the SCE decreases, the land surface albedo (a_s) decreases. This results in greater absorption of solar radiation in the Earth–atmosphere system, which reinforces further warming at the surface (Roesch and Roeckner 2006) and subsequent anomalies in the energy budget at the TOA (Flanner et al. 2011). Thus, the snow cover variation over the NH has gathered widespread attention because it is thought to be one of the most important factors in ongoing and projected climate changes. For example, McCabe and Wolock (2009) proved that long-term variability in the NH SCE in March has decreased substantially since 1970. Brown and Robinson (2011) used observational evidence to determine that the NH spring SCE has undergone significant reductions during the past ~90 years. Derksen and Brown (2012) further confirmed that spring SCE reductions in 2008–2012 have even exceeded the Coupled Model Inter-comparison Project Phase 5 (CMIP5) model projections.

In conjunction with the substantial attention on NH snow cover variation, methods have been developed to quantify the impacts of snow cover change. Flanner et al. (2011) assessed the influence of the changes in the NH cryosphere, including snow cover and sea ice, on the Earth’s radiation budget at the TOA. They concluded that the changes in cyrospheric (snow and ice) radiative forcing (C_rRF) from 1979 to 2008 was 0.45 Wm⁻² and changes in snow cover contributing about half of the total C_rRF over the NH. Chen et al. (2015) quantified the radiative forcing induced by snow cover phenology changes from 2001 to 2014 and found that the snow radiative forcing in the melting season accounts for 51 % of the total shortwave flux anomalies at the TOA in the corresponding period. Singh et al. (2015) estimated the global annual-mean land C_rRF during 2001–2013 with different kernels and different algorithmic determinations of snow presence and surface albedo and found that the global annual-mean land C_rRF during 2001–2013 was −2.6 Wm⁻². Perket et al. (2014) diagnosed current and projected shortwave cryosphere radiative effects and found that the global average NH cryosphere forcing from snow and sea ice is presently −4.0 Wm⁻², but the Earth will absorbs 1.6 Wm⁻² more insolation from snow and sea ice loss by 2099 according to Representative Concentration Pathways (RCP) 8.5 simulations.

Quantifying the impacts of snow cover changes is essential for meteorological, hydrological, ecological, and societal implications. However, unlike large-scale snow cover variation and its effects on the Earth’s climate system, the impacts of snow cover changes on the local scale have not been well studied, particularly in mainland China. As one of the most important components of the NH seasonal snow covered landmass, with the highest population and second-largest economy, snow cover condition in mainland China has drawn wide attention over the past decades (Qin et al. 2006; Shi et al. 2011; Shi and Wang 1980; Yang et al. 2008; Zhang et al. 2010). However, no consistent trend over the entire country has been identified in previous studies, and few have estimated the subsequent snow-induced radiative forcing. For example, Qin et al. (2006) reported that the snow cover during 1951–1997 showed a small increase trend over western China. Yang et al. (2008) found a deepening maximum snow depth at a rate of 1.15 cm decade⁻¹ during 1959–2003 in the Tianshan Mountains with a slight increase in the number of snow cover days. Zhang et al. (2010) monitored snow cover conditions in Liaoning Province, north-eastern China, but focused only on the winter months of November–April 2006–2008. Shi et al. (2011) demonstrated a decreasing trending in snow cover days and snow cover depth over China during 1990–2010 and predicted a further decreasing trend of snow cover days and snow depth until 2100. However, no consistent trend of snow cover can be drawn from these studies and few of these studies have discussed or analyzed the subsequent impacts of snow cover variation on the Earth’s climate system.

Under this background, the aim of the present study is to investigate the spatiotemporal variation in snow cover over mainland China and to quantify the associated snow-induced radiative forcing from 1982 to 2013 by using satellite observations with high spatial and temporal resolutions. Section 2 gives a brief review of the datasets and methods used in this study. The results are presented in Section 3, and a discussion and conclusions are given in Section 4.

2.1 Study area

In China, most of the snow is distributed in the mountains regions of the Tibetan Plateau (TP), Xinjiang, and northeastern China. In the distributions of domains and ecoregion maps (Xie et al. 2012), the snow-covered area covers four ecoregions including the northwestern arid (NWA) zone, the TP zone, the Inner Mongolia Plateau (IMP) zone, and the northeastern humid and semi-humid (NEH) zone. The four ecoregions and an elevation map of mainland China are shown in Fig. 1. To compare snow cover among different years, this study is confined to stable snow-covered regions (grid cells) where snow has covered the ground for at least 60 days in 80 % of the years during 1982–2013.

a The four ecoregions of mainland China. b The elevation map of mainland China

Full size image

2.2 Datasets

2.2.1 Snow cover datasets

To investigate the snow cover conditions at the highest achievable spatial and temporal resolutions, satellite-derived snow cover observations from both the National Oceanic and Atmospheric Administration/National Climatic Data Center (NOAA/NCDC) Climate Data Record (CDR) of the Northern Hemisphere Snow Cover Extent (NHSCE) dataset (Robinson et al. 1993) and the Moderate Resolution Imaging Spectroradiometer Satellite (MODIS) monthly snow cover fraction in Climate Modeling Grid (MOD10CM) products (Hall et al. 1995) were used in this study.

To compare the SCE information derived from NHSCE and MOD10CM consistently and to identify the snow cover variation in each pixel, the snow cover fraction (SCF) was used as the snow cover index in this study. The NHSCE dataset from 1982 to 2013 (Robinson et al. 1993) in a 25 km spatial resolution of roughly 0.25° was used to build the long-term series of SCF over China. The MODIS SCF products during 2001–2013 developed by Hall et al. (1995) in 0.05° spatial resolution were employed to generate a SCF map with finer spatial resolution. These products were used as the bench value of SCF series over mainland China because MODIS provides more consistent, objective snow estimates derived from high-resolution visible satellite data (Frei et al. 2012) compared with the NHSCE snow datasets used in this study. For the binary weekly NHSCE dataset, we reproduced monthly averaged and yearly averaged mean fractional SCE values as the SCF in each pixel, in which the SCF was defined as the total number of days with snow observations in a given period divided by the number of days in the period. For the MOD10CM monthly mean fractional SCE dataset, we used the fractional SCE value as the SCF value directly; SCF is the monthly mean fractional snow cover in a given pixel generated by using all of the days of a month.

2.2.2 Land surface albedo dataset

To quantify the land surface albedo contrast (∆a_s) applied to estimate the snow radiative forcing (RF_snow), which in this study is defined as the instantaneous perturbation to the Earth’s shortwave anomaly induced by snow cover changes, both the Global Land Surface Satellite (GLASS) a_s product (Liang et al. 2013) and the MODIS-derived global gap-filled snow-free a_s (MCD43GF) datasets (Moody et al. 2008) were used in this study. GLASS a_s is a real-time series of albedos, and the MCD43GF a_s is a snow-free albedo. Thus, we obtained the yearly and monthly Δa_s time series by subtracting MCD43GF snow-free a_s from GLASS a_s.

The GLASS a_s dataset provides a long-term series of high-quality, gap-filled global a_s maps from 1981 to 2013. GLASS a_s was produced from both Advanced very-high-resolution radiometer (AVHRR) and MODIS data (Liang et al. 2013) and has been used to quantify the radiative forcing of snow melting over Greenland (He et al. 2013) with high quality and fine spatial resolution. The spatial complete snow-free MCD43GF a_s is produced by using an ecosystem-reliant temporal interpolation technique that retrieves missing data with 3–8 % error (Moody et al. 2008).

To obtain accurate results of RF_snow in this study, we calculated both the white-sky and black-sky ∆a_s, respectively. In the reproduced processes, the selected quality control flags with “00” and “01” of GLASS a_s, indicating uncertainties of <5 and 10 %, were used to generate the monthly mean value of GLASS a_s. Quality assurance values between “0” and “1” of MCD43GF snow-free a_s, representing majority high-quality full inversion values, were used to generate the climatology value of monthly MCD43GF snow-free a_s.

2.2.3 Radiative kernels dataset

The albedo radiative kernel expressed as the TOA net shortwave anomalies associated with a 1 % change in the a_s estimated by 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) were used to quantify the RF_snow in this study. Compared with the analytical models, the radiative kernel methods are more accurate because they are more effective than the analytical models for capturing the functional dependence of planetary albedo (a_p) on a_s (Qu and Hall 2014). In this study, both the all-sky and clear-sky albedo radiative kernels were used to calculate the RF_snow in corresponding conditions.

2.2.4 Temperature and precipitation datasets

The surface air temperature dataset gridded at 0.25° horizontal resolution derived from the European Center for Medium Range Weather Forecasts Reanalysis (ERA-Interim) (Dee et al. 2011) and precipitation data set at 0.50° spatial resolution derived from the Climatic Research Unit (CRU) Time-Series (TS) Version 3.22 (Harris et al. 2013) during 1982–2013 were employed to attribute the snow cover variation over China. We did not use the precipitation dataset from ERA-Interim in our analysis mainly because its precipitation estimates are higher than those of the University of Washington, Asian Precipitation–Highly-Resolved Observational Data Integration Towards Evaluation of (APHRODITE) Water Resources, and Tropical Rainfall Measuring Mission Multi-satellite Precipitation Analysis datasets for all seasons and all basins in the TP zone (Tong et al. 2014).

To explore the contributions of the accumulation and melting season temperatures and precipitation to the snow cover changes over mainland China, we separated the entire snow season into an accumulation season and a melting season. On the basis of the seasonal cycle of snow cover (Fig. 2c), we defined the hydrological year of mainland China from September of 1 year (t − 1) to August of the following in year (t). Moreover, we defined the accumulation season of snow cover as September of 1 year (t − 1) to next February of the following year (t) and melting season as March to August in 1 year (t).

13-year averaged annual mean SCF over mainland China derived from a NHSCE and b MOD10CM datasets from 2001 to 2013. c Monthly mean SCF series over mainland China derived from NHSCE (blue line) and MODIS (red line) datasets from 2001 to 2013 and intergrated SCF series (green line) during 2010–2013. d Linear regression model in accumulation season and melting season between NHSCE and MODIS SCF series from 2001 to 2010. Y_a and Y_m in d represent the linear regression model in accumulation season and melting season, respectively. Correlations in d are significant at the 0.05 level

Full size image

Before calculation and analysis, all of the datasets used in this study were regridded at a spatial resolution of 0.05° with geographic latitude/longitude projection by using the resampling method of “average” or “cubicspline” with the help of gdalwarp (http://www.gdal.org/gdalwarp.html). For datasets with spatial resolution greater than 0.05°, we used “average” in the resampling process, which computed the average of all non-NODATA contributing pixels in the domain of our study. For datasets with spatial resolution lower than 0.05°, we used “cubicspline” in the resampling process.

2.3 Methods

2.3.1 Intergrated SCF series

The spatial distribution of annual mean SCF maps derived from NHSCE and MODIS were similar during 2001–2013 (Fig. 2a and 2b). However, compared with the MODIS SCF maps (Fig. 2b), the SCF values in the southern TP zone and NEH zones were overestimated in NHSCE SCF maps (Fig. 2a). The overestimation in these maps was caused mainly by low spatial resolution as well as the definition of snow cover in the NHSCE snow cover products, in which snow was marked as 1-byte values of 0 or 1, with <50 or ≥50 % probability of occurrence, respectively. The definition of snow cover in the NHSCE dataset may result in overestimated (underestimated) SCF values in regions of heavy (low) snow distribution because NHSCE can detect only pixels with ≥50 % probability of snow cover occurrence; the scattered distributed snow cannot be effectively identified.

To overcome the issues in the NHSCE SCF maps displayed in Fig. 2a and to build a long-term SCF series over mainland China, the NHSCE SCF maps were adjusted to fit the MODIS SCF maps by using a linear regression model. Both the NHSCE and MODIS datasets cover the period 2001–2013. Thus, the 2001–2010 monthly SCF series were used to build the linear regression model in each month, and the 2011–2013 SCF series were used to verify the reliability of linear regression model in each month.

The NHSCE monthly mean SCF series showed higher values in the peaks compared with the MODIS monthly mean SCF series in the intra-annual cycle of snow cover fluctuation during the hydrological years of 2001–2013 (Fig. 2c). However, except for some extremely high values, most of the plots generated from the MODIS SCF series can be well represented by the NHSCE SCF series by using the linear regression model, both in accumulation season and in melting season (Fig. 2d). The squared correlation coefficient (R²) between these two SCF series in accumulation season and melting season over the entire mainland China were 0.79 and 0.72, which are significant at the 0.05 level. By using the aforementioned linear regression model, we obtained the intergrated SCF series from 2011 to 2013 (Fig. 2c, green line). We found that the intergrated SCF series was highly similar to MODIS SCF series. The root–mean–square error (RMSE) between the intergrated SCF series and MODIS SCF series was 3.38 % during hydrological years of 2011–2013. Thus, it is reasonable to generate a long-term, intergrated monthly SCF series based on the NHSCE SCF series by using this linear regression model. By using this approach, we calculated the linear regression model in each grid cell and derived the intergrated annual and monthly SCF series.

2.3.2 Snow radiative forcing calculation

The radiative kernels method is widely used in the estimation of radiative forcing induced by a_s changes. In our study, the radiative kernels method was employed to quantify the RF_snow over mainland China, and we defined the RF_snow as the instantaneous perturbation to Earth’s TOA shortwave flux induced by the presence of snow cover. The time (t)-dependent RF_snow, within the study area R of area A composed of grid cells r can be represented as (Chen et al. 2015; Flanner et al. 2011)

Here, S is the spatial range (grid cells) of snow cover over the study area, ∂a_s/∂S is the dependence of a_s variation on the snow cover changes, and ∂F/∂a_s is the response of TOA flux variation to the a_s changes. We made an assumption following Flanner et al. (2011) such that monthly and spatially varying ∂a_s/∂S and ∂F/∂a_s were consistent with the snow cover fraction S and surface albedo a_s, respectively. Then, ∂a_s/∂S can be replaced with the mean ∆a_s observations that are specific to the applied snow cover data, and ∂F/∂a_s can be obtained from the albedo radiative kernels. The monthly gridded intergrated SCF, black-sky and white-sky climatologies of ∆a_s, and two albedo radiative kernels, CAM3 and AM2, in all-sky and clear-sky conditions were employed to calculate the RF_snow in this study.

2.3.3 Sensitivity analysis

SCF is largely determined by the snow season temperature and precipitation, including the accumulation season temperature (T_a) and precipitation (P_a) as well as the melting season temperature (T_m) and precipitation (P_m). To attribute the snow cover changes over mainland China from 1982 to 2013, we calculated the sensitivity of SCF to T_a, P_a, T_m, and P_m separately.

For each grid cell, we regressed the annual mean SCF as a dependent variable with T_a, P_a, T_m, and P_m as independent variables.

This formula was employed by Peng et al. (2013) to identify the sensitivity of the snow end date to the temperature and snow depth.

3.1 Snow cover variation over mainland China

The spatial distribution of the 32-year annual mean SCF and 32-year changes during 1982–2013 are displayed in Fig. 3a and 3b. The changes are expressed as linear slope multiplied by the time interval in this study. The 32-year annual mean SCF values were 7.43 ± 0.46 % over mainland China, with 16.86 ± 2.49 %, 10.85 ± 0.53 %, 8.35 ± 1.26 %, and 7.97 ± 0.97 % in NEH, NWA, IMP, and TP zones, respectively. The high SCF values were distributed mainly in the northern NWA, eastern NEH, and southern TP regions (Fig. 3a). The 32-year changes in SCF displayed large differences among the four zones, in which the SCF increased 3.73 %, 1.94 %, and 1.39 % in the NEH, TP, and IMP zones, respectively, but decreased 0.85 % in the NWA zone (Fig. 3b). In these regions, the maximum positive change of 31.08 % occurred in the southeast TP, and the maximum negative change of 27.49 % was distributed in northern Xinjiang and in the western margin region of the TP zone.

32-year a annual mean SCF, b changes, c inter-annual variability in annual mean SCF and linear trend in SCF, and d intra-annual variability in monthly mean SCF and 32-year changes of SCF in each month during 1982–2013. The black dots in b and solid bars in d indicate changes are significant at the 0.05 level, in which the orange (green) color indicates an increase (decrease) in SCF in a given month from 1982 to 2013. The error bar in d is expressed as the standard deviation of SCF in a given month during 1982–2013

Full size image

The 1982–2013 annual mean and monthly mean SCF over the entire mainland China calculated by the intergrated SCF maps are shown in Fig. 3c and 3d. The 32-year SCF in mainland China displayed an overall increasing trend of 0.29 % decade⁻¹ from 1982 to 2013 (Fig. 3c), which is significant at the 0.05 level. Detailed analysis of inter-annual variability of SCF in mainland China determined by the Mann–Kendall test indicated a decreasing trend of −0.12 % decade⁻¹ before 1996 and an increasing trend of 0.37 % decade⁻¹ during 1997–2013 with peaks occurring in 2000 and 2010 and valleys presented in 2002 and 2008.

The intra-annual cycle of SCF presented an increase (decrease) in SCF in January (June) of 4.13 % (−1.23 %) from 1982 to 2013, which is significant at the 0.05 level (Fig. 3d). Moreover, the 32-year changes in the intra-annual cycle displayed an increase in SCF from October to the following May and a decrease in SCF from June to September, which indicates that the intra-annual amplitude of SCF increased during 1982–2013 and that more snow cover was released as snowmelt runoff in spring.

3.2 Attribution of SCF changes

The variation in snow cover is closely related to changes in temperature and precipitation, which is driven by large-scale atmosphere circulation. A seasonal increase in the mean temperature of 3 to 5 K reduces the snow cover duration by more than 1 month on average (Keller et al. 2005). To attribute the SCF variation over mainland China, the 32-year changes in T_a, T_m, P_a, and P_m were calculated from the linear regression trend, as shown in Fig. 4a to 4d. The sensitivity of SCF to changes in T_a, T_m, P_a, and P_m calculated from Eq. (2) are shown in Fig. 4e to 4h.

32-year changes in a T_a, b T_m, c P_a, and d P_m over mainland China during 1982–2013. e Sensitivity of SCF to T_a when the impacts of T_m, P_a, and P_m to SCF are removed. f Sensitivity of SCF to T_m when the impacts of T_a, P_a, and P_m to SCF are removed. g Sensitivity of SCF to P_a when the impacts of T_a, T_m, and P_m to SCF are removed. h Sensitivity of SCF to P_m when the impacts of T_a, T_m, and P_a to SCF are removed. The black dots in a to d indicate changes are significant at the 0.05 level. Those in e to h indicate linear correlation coefficients between SCF and T_a in e, SCF and P_a in f, SCF and T_m in g, and SCF and P_m in h are significant at the 0.05 level

Full size image

As displayed in Fig. 4a, most of regions over mainland China experienced a cooler T_a from 1982 to 2013, which is beneficial for snow accumulation. The decreased T_a is essentially consistent with the long-term tendency of large-scale cooling trends of land surface temperature in winters over mid-latitudes observed since about the 1990s (Cohen et al. 2012, 2014). Changes in the Arctic system were attributed to a T_a decline in mid-latitudes through changes in storm tracks, the Jet Stream, and planetary waves and their associated energy propagation (Cohen et al. 2014). Moreover, the recent Pacific Ocean cooling effect also contributed to the declined T_a in mid-latitudes because the tropical cooling effect on the extra-tropics is most pronounced in winter (Kosaka and Xie 2013). In contrast to changes in T_a, most of the regions over mainland China experienced a warmer T_m from 1982 to 2013 (Fig. 4b), especially the NWA and NEH zones.

Accompanied by the decline in T_a in mainland China, the P_a increased in most of northern China and the TP regions from 1982 to 2013 (Fig. 4c), which is beneficial to snow accumulation and may have resulted in deeper snow depth and subsequently a longer snow season. This pattern coincides with the previous findings of large-scale cold snaps, heavy snowfall, and glacier events at middle latitudes (Cohen et al. 2010, 2012; Yao et al. 2012). The increased precipitation at such latitudes is attributed to afforestation (Swann et al. 2012) as well as strengthened Westerlies (Mölg et al. 2013; Yao et al. 2012). The cooling trend of T_a accompanied by excessive precipitation from 1982 to 2013 result in an increased SCF in the TP regions, which is similar with the difference of glacier conditions distributed in the TP and surroundings (Yao et al. 2012). The changes in P_m (Fig. 4d) differed significantly from the P_a anomalies during the same time span, with increased P_m in the TP and NWA zones and decreased P_m in the NEH regions. The positive anomalies of P_m in the TP and NWA zones agree with previous study results that show a change in the climate in northwest China from dry-arid to warm-humid since 1987 (Shen et al. 2013).

The sensitivity analysis in Eq. (2) displayed a negative sensitivity of SCF to T_a, at −0.75 %°C⁻¹, and a positive sensitivity to P_a, at 0.34 % cm⁻¹, over the entire area of mainland China from 1982 to 2013 (Fig. 4e and 4g). Compared with the sensitivity of SCF to temperature and precipitation anomalies in the accumulation season, the magnitudes of those in melting season were significantly smaller (Fig. 4f and 4h), at −0.12 %°C⁻¹ and 0.07 % cm⁻¹. However, the correlation analysis between SCF and T_a as well as P_a demonstrated that the SCF was largely determined by T_a anomalies (R² = 0.79), as P_a was not statistically correlated with the SCF (R² = 0.18) at the 0.05 significance level over mainland China except for the NEH zone.

3.3 Snow radiative forcing

To reduce the uncertainty in the RF_snow calculation, we averaged two RF_snow results using CAM3 and AM2 albedo radiative kernels seperately to obtain the mean RF_snow values. The spatial distributions of the 32-year annual mean clear-sky RF_snow and those changes during 1982–2013 are shown in Fig. 5a and 5b. The inter-annual and intra-annual variability of clear-sky RF_snow and the associated changes over the entire mainland China are shown in Fig. 5c and 5d.

32-year a annual-mean RF_snow, b changes, c inter-annual variability of annual-averaged RF_snow and linear trend in RF_snow and d intra-annual variability of monthly averaged RF_snow and 32-year changes of RF_snow in each month over mainland China during 1982–2013. The annual-mean values in a and c were averaged from two RF_snow results using CAM3 and AM2 kernels. The error bars in c and d are expressed as the range derived from the two RF_snow results in clear-sky condition. The solid green bars in d indicate that the RF_snow was enhanced, and the changes are significant at the 0.05 level in a given month over mainland China during 1982–2013

Full size image

The 32-year annual mean value of clear-sky RF_snow over the entire mainland China during 1982–2013 was −4.33 ± 0.29 Wm⁻² with −8.46 ± 1.19 Wm⁻², −5.85 ± 0.49 Wm⁻², −5.67 ± 0.71 Wm⁻² and −5.79 ± 1.01 Wm⁻² recorded in the NEH, NWA, IMP, and TP zones, respectively. The large values of clear-sky RF_snow were distributed mainly in the southeastern TP, northern NWA, and northeastern NEH zones because of the high annual mean SCF values in these regions. In contrast, the changes in clear-sky RF_snow were distributed in the southeastern and southwestern TP zones, the former of which showed a significant RF_snow change of −8.26 Wm⁻² during 1982–2013. This result indicates that the cooling effects caused by snow cover have been enhanced since 1982. In contrast, obviously weakened cooling effects were distributed in the northeastern NEH zones and southwestern TP which coincides with the footprint of decreased SCF in this region (Fig. 3b). Induced by the weaker Indian monsoon (Yao et al. 2012), the winter precipitation decreased in the southwestern TP zone, which resulted in a decline in SCF followed by a lower ∆a_s and reduced RF_snow in this region.

Compared with the 32-year annual mean SCF anomaly displayed in Fig. 3b, it is evident that the fluctuation in RF_snow is generally opposite the changes in SCF, indicating that the increases in SCF relate to a greater ∆a_s and additional solar radiation reflected to space, which resulted in a strengthened cooling effect of snow cover during 1982–2013. The 32-year change in clear-sky RF_snow over mainland China during 1982–2013 was −0.65 Wm⁻² with −1.44 Wm⁻², −1.15 Wm⁻², −0.65 Wm⁻², and −0.58 Wm⁻² recorded in the IMP, TP, NWA, and NEH zones, respectively. The peak of annual mean clear-sky RF_snow, at −3.47 (±0.09) Wm⁻² occurred in 2007, whereas the valley, at −5.71 (±0.14) Wm⁻² occurred in 2005.

The intra-annual cycle of clear-sky RF_snow and 32-year linear changes of clear-sky RF_snow in each month are shown in Fig. 5d. The clear-sky RF_snow peaked broadly at approximately −11.19 (±0.17) Wm⁻² in February–March because of the high value of solar radiation and SCF in that month. The clear-sky RF_snow valleys were broad at approximately −0.09 (±0.001) Wm⁻² during July–September mainly because of the low SCF as well as low albedo radiative kernels owing to higher cloud optical thickness over the NH compared with that in other months. During October–January, the clear-sky RF_snow enlarged owing to the increase in SCF, although the insolation was still low in these months. Except for July and August, all of the months showed negative changes, particularly in March in early spring and November in early winter.

The all-sky RF_snow showed inter-annual, intra-annual, and spatial pattern variabilities similar to those in clear-sky conditions (the result is not displayed in the figures). The 32-year annual mean value of all-sky RF_snow over mainland China during 1982–2013 was −3.01 ± 0.18 Wm⁻². The peak and valley values of annual mean all-sky RF_snow over mainland China were −2.29 (±0.01) Wm⁻² in 2007 and −3.76 (±0.02) Wm⁻² in 2005, respectively, consisting a 30 % difference in clear-sky values. Compared with the estimates of annual mean RF_snow in clear-sky and all-sky conditions with white-sky albedo, those using black-sky albedo showed a slightly larger value of 4.32 %. This this difference has also been reported by Flanner et al. (2011) and Qu and Hall (2014).

On the basis of satellite observations at high spatial and temporal resolutions, this study has described the distribution of snow cover changes over mainland China and their attribution and has quantified the associated RF_snow from 1982 to 2013. These results can help to provide a better understanding of the spatial and temporal changes of snow cover over mainland China in the context of climate change to more clearly evaluate its contribution to the Earth’s climate system.

In contrast to the widely known snow cover reduction over the NH (Brown and Robinson 2011; Derksen and Brown 2012; McCabe and Wolock 2009), the snow cover in mainland China has increased from 1982 to 2013 at a rate of 0.29 % decade⁻¹, which is significant at the 0.05 level. The increased SCF over mainland China during 1982–2013 agrees with previous findings of large scale cold snaps, heavy snowfall, and glacier events across the United States, Europe, and East Asia and the observed long-term tendency of large-scale cooling trends of land surface temperature in winter over mid-latitudes. In addition, our results are consistent with previous reports of increased SCE over western China (Qin et al. 2006) and deepened maximum snow depth in the Tianshan Mountains (Yang et al. 2008). However, they disagree with the simulation results of decreasing trends of snow cover days and snow cover depth over mainland China during 1990–2010 reported by Shi et al. (2011). Accompanied by the snow cover variation, the RF_snow enhanced 0.21 Wm⁻² decade⁻¹ during the corresponding period, which is opposite the weakened cooling effects of NH snow and sea ice (Flanner et al. 2011).

Unlike previous research in snow cover, this study used both NHSCE and MODIS snow cover datasets in the investigation of snow cover variation over mainland China, which ensured the reliability of snow cover quantification. In addition, this study used GLASS and MCD43GF a_s products in the calculation of RF_snow, which greatly reduced the uncertainty in RF_snow estimation caused by a_s biases that exists in previous studies. The snow-covered a_s used by Flanner et al. (2011) was used in decreasing priority with the monthly-resolved MODIS a_s product, then its annual-mean; the AVHRR Polar Pathfinder (APP)-x data; and MODIS land-class means. Thus, the spatially complete GLASS a_s is more effective for RF_snow estimation. However, the unsolved spectral variation in the radiative kernel technique and effects of vegetation in the RF_snow estimation were not included in this paper. Snow induces greater ∆a_s in the visible portion of the spectrum than that in the near-infrared, whereas the radiative kernels are generated with spectrally uniform surface albedo perturbations but with spectrally varying surface fluxes. However, sensitivity analysis by Flanner et al. (2011) found that the unsolved spectral variation in the radiative kernel technique may result in a 5 % difference when comparing C_rRF estimates derived from solar broadband quantities with those derived from the same products partitioned into visible and near-infrared components. In addition, changes in vegetation type may effect snow-free a_s, and result in the corresponding uncertainties in RF_snow calculation (Flanner et al. 2011; Harvey 1988). Compared with forest land, grassland may result in larger ∆a_s when snow occurs (Harvey 1988). However, study by Harvey (1988) proved that the radiative forcing change caused by land snow cover changes is about three times larger in the absence of vegetational masking over the NH. Thus, the excluded unresolved spectral variation in the radiative kernel technique and the effects of vegetation in RF_snow estimation did not influence the main conclusions of this study.

The global cold-season warming has enhanced in semi-arid regions from 1901 to 2009 (Huang et al. 2012). However, cold-season cooling in the arid and semi-arid NWA and TP zones of China increased from 1982 to 2013 under this background. This result is opposite that for long-term warm winters over the NH. The increased SCF over mainland China benefits regions in which the water supply depends on snowmelt, particularly in arid and semi-arid regions. However risks are present because floods may be generated by the melting of accumulated snow in spring (Arnell and Gosling 2014). Thus, it is also important to monitor and predicate the peaks and volumes of snow-induced runoff over mainland China, particularly in northwest arid and semi-arid mountain regions. In addition, global climate change is expected to result in greater variation in snow cover and subsequent impacts on land surface hydrology and vegetation production (Paudel and Andersen 2012). Thus, the contribution of snow cover in areas of vegetation change and the vegetation effects on RF_snow estimation are also worthy topics for future research.