Quantifying Snow Albedo Radiative Forcing and Its Feedback during 2003–2016

Abstract

Snow albedo feedback is one of the most crucial feedback processes that control equilibrium climate sensitivity, which is a central parameter for better prediction of future climate change. However, persistent large discrepancies and uncertainties are found in snow albedo feedback estimations. Remotely sensed snow cover products, atmospheric reanalysis data and radiative kernel data are used in this study to quantify snow albedo radiative forcing and its feedback on both hemispheric and global scales during 2003–2016. The strongest snow albedo radiative forcing is located north of 30°N, apart from Antarctica. In general, it has large monthly variation and peaks in spring. Snow albedo feedback is estimated to be 0.18 ± 0.08 W∙m⁻²∙°C⁻¹ and 0.04 ± 0.02 W∙m⁻²∙°C⁻¹ on hemispheric and global scales, respectively. Compared to previous studies, this paper focuses specifically on quantifying snow albedo feedback and demonstrates three improvements: (1) used high spatial and temporal resolution satellite-based snow cover data to determine the areas of snow albedo radiative forcing and feedback; (2) provided detailed information for model parameterization by using the results from (1), together with accurate description of snow cover change and constrained snow albedo and snow-free albedo data; and (3) effectively reduced the uncertainty of snow albedo feedback and increased its confidence level through the block bootstrap test. Our results of snow albedo feedback agreed well with other partially observation-based studies and indicate that the 25 Coupled Model Intercomparison Project Phase 5 (CMIP5) models might have overestimated the snow albedo feedback, largely due to the overestimation of surface albedo change between snow-covered and snow-free surface in these models.

1. Introduction

The globally averaged surface temperature has increased by 0.85 °C over the period of 1880–2012 [1], with particularly strong warming signals appearing at high northern latitudes [2,3,4]. This is known as the Arctic amplification [5,6,7]. One of the main contributions to the amplified warming is perhaps the positive surface albedo feedback [8,9,10], primarily snow and ice albedo feedback. The ice, specifically, refers to bare ice, melting ice (mainly includes sea ice, glaciers and ice caps), snow-covered ice, open water, etc. [1]. In the warming climate, the decreasing snow cover (ice) extent and snow (ice) depth are leading to a less reflective planet that absorbs more solar radiation, and thus warming the earth further [11,12].

Climate feedback variables are valuable indicators of climate change due to their sensitivity to temperature. Specifically, a central task of climate change research is to quantify the Equilibrium Climate Sensitivity (ECS, [13,14,15]), which refers to the equilibrium change in annual mean surface temperature with a doubling of the atmospheric equivalent carbon dioxide concentration [16,17,18]. However, no agreement has been reached on the magnitude of the ECS [19], and the large uncertainty in the ECS is primarily attributed to the inaccurate estimation of individual feedbacks [20,21,22].

Models are useful tools for climate feedback study for their capability of long-term simulation. Energy balance models were first used mainly for the recognition of feedback mechanisms and dynamic processes [23,24,25]. Later on, General Circulation Models (GCMs) were widely used to investigate whether variables act as positive or negative feedbacks and the magnitude of each feedback [26,27,28]. In recent development, the Atmosphere–Ocean General Circulation Models (AOGCMs) make it possible for more detailed and comprehensive estimation of individual feedbacks [17,29,30]. However, feedback estimations from model simulations are still associated with large uncertainties, i.e., fivefold intermodel difference on surface albedo feedback was reported in the Fourth Assessment Report (AR4) of Intergovernmental Panel on Climate Change (IPCC) models [31]. Such uncertainty remains in IPCC AR5 [1,32].

Observation-based research has also made its contribution to feedback assessment, with more realistic representation of the climate from observations. For example, Flanner et al. found that the Coupled Model Intercomparison Project Phase 3 (CMIP3) models underestimated the snow and ice albedo feedback in Northern Hemisphere, as compared with results based on satellite observations of the Extended Advanced Very High Resolution Radiometer Polar Pathfinder (APP) product and Moderate Resolution Imaging Spectroradiometer (MODIS) data [33]. Fletcher et al. showed that the modeled average snow albedo feedback in the CMIP3 models was slightly larger than satellite observations [34]. Dessler reported an agreement on both global average and spatial pattern between model results and individual feedbacks that were calculated based on reanalysis data of ERA-Interim [35] and Modern-Era Retrospective Analysis for Research and Applications (MERRA, [36]) [37]. However, problems also exist in observation-based results. For example, the coarse resolution of observational data would either temporally or spatially smooth the feedbacks, especially in snowmelt seasons [38,39].

Due to the fact that large differences and uncertainties remain in snow albedo feedback assessments, an effective way to improve its accuracy is to use observation-based results to constrain model assessments [1]. Nonetheless, most snow albedo feedback studies were based on model simulations, and only few on observations [40,41,42]. Moreover, most of these studies calculate the combining effects of snow albedo feedback and ice albedo feedback, namely surface albedo feedback [13,37,43]. As a result, it is impossible to separate the contribution of snow or ice albedo feedback, as well as their uncertainties.

Being motivated by these scientific findings and limitations, this study specifically focuses on the quantification of snow albedo feedback by remotely sensed snow cover products of MODIS (MOD10C1 and MYD10C1), atmospheric reanalysis data and radiative kernel data. The purpose is to examine the source of the differences between our result and partially observation-based results, as well as Coupled Model Intercomparison Project Phase 5 (CMIP5) model-based results, and to provide precise information of snow cover change and snow albedo radiative forcing for model parameterization.

Data and method are described in Section 2. Spatial and temporal variability of snow cover, snow albedo radiative forcing and its feedback are calculated and analyzed in Section 3. Comparison and discrepancy analysis with both partially observational-based and model-based studies, strengths and uncertainties of this work are discussed in Section 4. Conclusion follows in Section 5.

2. Materials and Methods

2.1. Data

2.1.1. Remote Sensing Data

MODIS is a key imaging instrument onboard two complementary satellites, Terra and Aqua, and provides two series of datasets with different passing times each day. With relatively high spatial, temporal and spectral resolutions, MODIS offers timely monitoring for land, ocean, and atmosphere research [44]. In this paper, daily L3 Climate Modeling Grid (CMG) of fractional snow cover products from Terra (MOD10C1) and Aqua (MYD10C1) were selected.

The two products are with the same spatial resolution of 0.05°. MOD10C1 products were available since February 2000, and MYD10C1 products from July 2002. Both products are still being updated. Therefore, data duration of this study is chosen from 1 September 2002 to 31 August 2016. During the study period, if either image of Terra (MOD10C1) or Aqua (MYD10C1) were missing, the alternative image would represent the data of the day. If both images are missing (occurs only once during the study period), the mean value of the previous day and the latter day is considered as the missing day’s snow cover. We define every snow year from 1 September of the previous year to 31 August of the current year. Seasons are defined as autumn (September–October–November), winter (December–January–February), spring (March–April–May) and summer (June–July–August). Note that the spatial resolution of snow cover data were resampled to 0.25° when calculating snow albedo radiative forcing and snow albedo feedback, while it remained 0.05° in snow cover change analysis.

A general disadvantage of snow cover optical remote sensing products, e.g., MODIS snow cover products, is that clouds shadow the information of their underlying surface, usually resulting in an underestimate of snow cover. In addition, cloud changes in position and extent at the different passing time with different viewing angle of Terra and Aqua, thus the two sensors possibly detect different snow cover information. Therefore, daily combination method was used to get the potentially maximum snow cover against the block of cloud, which could be a more realistic representation of the actual snow cover amount [45]. The combination method combines snow cover information of Terra and Aqua at daily scale, and contains two circumstances. Firstly, if either image of Terra or Aqua in the same day considers the pixel to be snow, the pixel is considered as snow. Secondly, if both images of Terra and Aqua consider the pixel to be snow, but with different fractional snow cover value, the larger one is used as the fractional snow cover of the day.

2.1.2. Atmospheric Reanalysis Data

Atmospheric reanalysis data were extracted from ERA-Interim [35], which is the latest product of European Centre for Medium-Range Weather Forecasts (ECMWF). ERA-Interim is a global atmospheric reanalysis product that started in 1979, and is continuously updated in near real time [46].

Snow albedo and albedo data were extracted on global scale, with spatial and temporal resolutions of 0.25° and 6 h, respectively. Specially, the albedo data in ERA-Interim refer to the background snow-free surface albedo. Monthly mean air temperature data at 2 m height were also acquired to estimate temperature change during the study period, with the same spatial resolution of 0.25°. All these datasets were obtained for the period 2003–2016. Snow albedo and snow-free albedo datasets (four per day) were separately averaged to daily resolution.

2.1.3. Radiative Kernel Data

Radiative kernel data [47,48] were obtained from Community Atmosphere Model version 3 (CAM3, [49]) of National Center for Atmospheric Research (NCAR) (Boulder, CO, USA). Surface albedo radiative kernel under all sky (K_α) was obtained, describing the response of net shortwave radiation at top of the atmosphere (TOA) to a 1% additive increase in surface albedo [47]. K_α is a function of latitude, longitude, and month of year. The spatial resolution of K_α is approximately 2.8° on average, and varies a little with latitude. The data were then bilinear interpolated to 0.25° for spatial consistency with other products.

Spatial distribution of monthly surface albedo radiative kernel is illustrated in Figure 1. Monthly variability of K_α is strongly affected by the movement of the sun. High values mostly appear in low latitudes where there is large incident solar radiation. K_α decreases as latitude increases, and diminishes to 0 in areas with polar darkness. Moreover, as clouds act to mask the effects of changes in the albedo of the underlying surface, K_α has the lower values in persistent cloud-covered regions, (e.g., mid-latitude storm track area), and relatively higher values in cloud-free areas (e.g., low-latitude desert area) [47].

2.2. Method

For the advantage of less computationally expensive and easier comparison with other results, radiative kernel method is widely used in snow albedo radiative forcing and feedback estimations [50,51,52].

Feedback estimations using radiative kernel method usually decompose feedback into two factors, radiative kernel and climate response pattern. The former depends only on the radiative algorithm and base climate. The latter describes the change in mean climatology of the feedback variable between the two climate states [47]. Thus, the magnitude of feedback is simply the product of the two factors.

However, on purpose of getting more process information of feedback, snow albedo radiative forcing was calculated first in our study, using Equation (1):

$$G_s(t,R),\{\begin{array}{ll}0,& ,S(t,r)=0\\ \frac{{{\displaystyle \int }}_R^{}\frac{\partial \alpha }{\partial S}(t,r)\times \frac{\partial H}{\partial \alpha }(t,r)\times dA(r)}{A(R)}& ,S(t,r)>0\end{array}$$

(1)

Here, G_s(t,R) is snow albedo radiative forcing that describes the instantaneous influence of snow cover on TOA energy budget, with unit of W∙m⁻² [33]. G_s(t,R) is a function of time t and study area R, and R has a total area of A and grid size r. S(t,r) represents snow cover at time t and pixel r. Being the prerequisite of snow albedo radiative forcing, snow cover of each pixel is checked before calculation. According to Equation (1), if pixel r is snow-free at time t, i.e., S(t,r) = 0, no snow albedo radiative forcing occurs, i.e., G_s(t,R) = 0. Only pixels with a fractional snow cover larger than 0 contribute to the snow albedo radiative forcing. The magnitude of G_s(t,R) depends on the albedo contrast between snow-covered and snow-free circumstances, as well as the radiative kernel data. Specifically, $\frac{\partial \alpha }{\partial S}(t,r)$ is surface albedo change induced by snow cover change, which can be simplified as albedo contrast between snow-covered and snow-free circumstances ([33,53]). $\frac{\partial H}{\partial \alpha }(t,r)$ is the radiative change against the coincident change of albedo, namely the surface albedo radiative kernel. Note that $\frac{\partial \alpha }{\partial S}(t,r)$ and S(t,r), $\frac{\partial H}{\partial \alpha }(t,r)$ and $\alpha$(t,r) are assumed independent from each other, and to avoid the inconvenience of negative numbers, the absolute value of K_α is used in the calculation.

Then, snow albedo feedback can be quantified by the amount of additional net shortwave radiation at TOA as surface albedo decreases in association with a 1 °C temperature increase [11]:

$$\mathsf{\lambda }=\frac{∆\mathrm{Gs}(t,R)}{∆T(t,R)}$$

(2)

where λ is the feedback parameter, here considered as snow albedo feedback, with unit of W·m⁻²·°C⁻¹. ∆G_s(t,R) and ∆T(t,R) are monthly anomalies of G_s(t,R) and T(t,R) during the study period, respectively, and both are functions of time t and study area R. In this study, t is set to be monthly, and when calculating the regional mean, R is defined as landmasses north of 30°N, namely the North Hemisphere Extratropical Land (NEL). Accordingly, G_s(t,R) and T(t,R) are the monthly snow albedo radiative forcing and monthly mean surface air temperature averaged over the NEL, respectively. Note that G_s(t,R) and T(t,R) were area weighted [34] during the calculation process.

The coefficient of least square fit of ∆G_s(t,R) and ∆T(t,R) represents the magnitude of snow albedo feedback [54]. Considering the relatively short study period, i.e., 14 years with 168 monthly anomalies, and the small magnitude of both ∆G_s(t,R) and ∆T(t,R), the result may be subjected to substantial uncertainty due to random variations of G_s(t,R) and T(t,R). Therefore, in order to get a precisely snow albedo feedback and confidence level, the block bootstrap test was used. Specifically, the method considers data of each year (12 monthly ∆G_s(t,R) and ∆T(t,R)) as a block, thus 14 blocks of data are contained in the original dataset. The process includes three steps:

Step 1: First, randomly pick one block of data from the 14 blocks, pick another from the 14 blocks of data (there is possibility that the same block is chosen again), pick a third block, etc. Repeat this process until all the 14 blocks are included in the newly picked dataset, in which, some blocks of data may appear more than once.

Step 2: Evaluate snow albedo feedback of the newly picked dataset based on least square fit.

Step 3: Repeat Step 1 and Step 2 for a large number of times, 10,000 times in our case. Estimate the mean snow albedo feedback and its confidence level according to the 10,000 snow albedo feedback results.

3. Results

3.1. Spatial and Temporal Variability of Snow Cover

Through the daily combination method, annual mean fractional snow cover from 2003 to 2016 increased 3.94% and 4.69% from the original MOD10C1 and MYD10C1 datasets, respectively. Being the original dataset of MOD10C1, MOD10A1 has an absolute overall accuracy of ~93% [55]. Nevertheless, accuracy decreases when cloud is taken into consideration. Due to the large uncertainties in accuracy assessment of MOD10C1 (MYD10C1), i.e., accuracy difference between MOD10A1 and MOD10C1, cloud variabilities, as well as the uncertainties in contribution of snow cover under cloud to snow albedo radiative forcing and snow albedo feedback, the contribution of the accuracy improvement through daily combination method is not discussed.

Global mean fractional snow cover during 2003–2016 (Figure 2) is calculated as the ratio between the sum of daily fractional snow cover and total days of the study period. With large spatial variability, snow cover mainly distributes in landmasses north of 30°N, excluding Antarctica. Greenland is the most densely snow-covered area over the NEL where most of the land is covered by snow and ice all year round. The Tibetan Plateau should also be mentioned for its abundant snow cover in such low latitude. Furthermore, the whole Antarctica continent is deliberately mapped as snow in MODIS snow cover products, thus it is not discussed in the paper. Apart from Antarctica, snow cover can be seen in limited areas of the Southern Hemisphere, i.e., the Andes Mountains, southeast corner of Australia, and a few parts of New Zealand.

As MODIS is an optical sensor, snow-covered areas inside of the Arctic Circle during polar night cannot be detected. This would certainly result in an underestimate of the total snow amount and maybe some discrepancies on trend analysis. However, as solar radiation is the precondition of snow albedo radiative forcing and feedback, areas in darkness receive no insolation, consequently exert no radiative forcing and feedback to the climate system. In other words, the “missing data” have little influence on radiative forcing and feedback assessment, thus the detected snow cover can be regarded as effective snow cover in our study.

Time series of monthly and annual mean fractional snow cover over the NEL during 2003–2016 are presented in Figure 3a. Here, the NEL rather than the globe was chosen for analysis, because snow cover over the NEL makes up for about 98.45% of the total snow amount during the study period, apart from Antarctica. Monthly variations are large in a sense that snow cover from winter and spring months accounts for over 75.28% of the total snow amount on average during 2003–2016, while summer months contributes only about 6.98%. Specifically, March has the largest fractional snow cover through the year. February, April and January follow and the four months have a much larger fractional snow cover value above the average. On the contrary, apart from summer months, September has the smallest snow cover amount, and its fractional snow cover amount is similar to that of June. Annual mean fractional snow cover over the NEL shows small interannual variability, while there is a slight, but insignificant (at the p = 0.05 level) decline in total. Months with large fractional snow cover (e.g., March, February, April, and January) are also the main contributors to interannual variability. Inconsistency tendencies of interannual variabilities are found in months, e.g., January and November experience an increasing trend of fractional snow cover in the recent three years while a decreasing trend are found in most of the rest months.

The climatological monthly mean fractional snow cover over the NEL during 2003–2016 is displayed in Figure 3b. The amount of the effective snow cover presents a single-peaked shape (peaks in March) through the year. Effective snow cover increases continuously from September to March, indicating a successive process of accumulation, during which winter months are the strongest and fastest accumulation period. Accordingly, snow cover decreases from March to July (the August value being slightly higher than the July value), and the fastest ablation occurs in spring from March to May.

3.2. Snow Albedo Radiative Forcing

Annual mean snow albedo radiative forcing during 2003–2016 is shown in Figure 4. Areas with strong snow albedo radiative forcing also appear over the NEL, excluding Antarctica. However, areas with large fractional snow cover do not always represent strong snow albedo radiative forcing, and the peak value appears in the Tibetan Plateau, rather than Greenland. This is partly because the Tibetan Plateau is located at the relatively low latitude with much more insolation than the higher latitudes, and snow there can persist until relatively late spring. In addition, areas in western part of the U.S., northern part of Canada, mid-high latitudes of Russia, and the Andes Mountains along the west coast of South America also exhibit strong snow albedo radiative forcing.

Climatological monthly mean snow albedo radiative forcing over the NEL during 2003–2016 is displayed in Figure 5a. Standard deviations are shown as whiskers, which represent interannual variability of each month during the study period. Autumn months from September to November exhibit small snow albedo radiative forcing. Because snow cover over the NEL starts to accumulate in autumn (Figure 3b), there is neither much snow cover nor strong insolation. Winter experiences the largest expansion of snow cover extent and the smallest insolation among the four seasons. Snow albedo radiative forcing is much larger than that in autumn, and it increases from December to February. Snow albedo radiative forcing is the largest in spring months and peaks in April, as a result of both large snow cover extent and strong solar radiation. Spring also experiences the largest decrease of snow cover extent throughout the year. This leaves summer with very little snow cover amount. As a result, in spite of its largest insolation, snow albedo radiative forcing in summer is much smaller than that in spring.

April exhibits the largest snow albedo radiative forcing throughout the year, and its interannual variability is displayed in Figure 5b. There is an overall insignificant (at p = 0.05 level) decline in April snow albedo radiative forcing over the NEL during 2003–2016. Large interannual variability with a continuous decline during the last 4 years are experienced during the study period. In general, April snow albedo radiative forcing decreased about 6.10% during 2003–2016.

3.3. Snow Albedo Feedback

Snow albedo feedback estimation based on the block bootstrap test is explained in detail with an example trial, according to the method introduced previously (Section 2.2). Specifically, following Step 1, it takes 27 times of random choice before all the 14 blocks of data are picked. Details of each pick are displayed in

Table 1. Every Pick of the example test.

Pick ID	1	2	3	4	5	6	7	8	9	10	11	12	13	14
Block ID	07	13	06	04	09	06	10	08	08	04	13	16	11	16
Pick ID	15	16	17	18	19	20	21	22	23	24	25	26	27
Block ID	03	05	12	03	14	10	14	03	10	14	07	11	15

. Specifically, Pick ID represents the count of each pick, while Block ID is the block of data of the corresponding year. During the experiment, it is not until the occurrence of Block 15 (the 27th pick) that all 14 blocks of the original dataset are included. Thus, 27 blocks with 324 data are contained in the newly generated dataset.

Following Step 2, snow albedo feedback can be quantified based on the newly picked dataset. Scatterplot of the 324 $\mathsf{\Delta }$G_s(t,R) and ∆T(t,R) data are shown in Figure 6a. Colors of the dots represent the frequency of their occurrence (corresponding to

Table 2. Frequency of each block of data based on the example test.

Block ID	03	04	05	06	07	08	09	10	11	12	13	14	15	16
Frequency	3	2	1	2	2	2	1	3	2	1	2	3	1	2

). The dash line is the least square fit of the 324 data, i.e., snow albedo feedback strength of this trial.

Following Step 3, the mean snow albedo feedback and the 95% confidence level can be estimated after repeating the experiment for 10,000 times. The histogram of the 10,000 snow albedo feedback is shown in Figure 6b. The value of the red line refers to the mean snow albedo feedback, in our case, 0.18 ± 0.08 W·m⁻²∙°C⁻¹ over the NEL during 2003–2016. This implies that as surface albedo decreases in association with 1 °C temperature increase, snow albedo feedback would cause an increase of 0.18 ± 0.08 W∙m⁻² in the net shortwave radiation at TOA, averaged over the NEL.

By rescaling the NEL result with product of two factors: the ratio of the NEL area to the global area and the ratio of annual mean surface air temperature change of the NEL to global mean change [56], global snow albedo feedback can be scaled as 0.04 ± 0.02 W∙m⁻²∙°C⁻¹. However, due to the limited sample amount, the result may be subjected to substantial uncertainty.

4. Discussion

4.1. Comparison with Partially Observation-Based Studies

There are extremely limited observation-based snow albedo feedback studies, so only three of the compatible ones (with similar study area of the North Hemisphere) are chosen for further comparison: snow and ice albedo radiative forcing estimation by Flanner et al. for 1979–2008 [33], both snow albedo radiative forcing and feedback estimations by Chen et al. for 1982–2013 [57] and snow albedo feedback estimation by Peng et al. for 1980–2006 [58]. It should be noted that, although feedback was also quantified in Flanner and coworkers’ study, the regional mean value of the combining snow and ice albedo feedback makes it incomparable with our snow albedo feedback results.

The magnitude of snow albedo feedback is examined first. By using a linear fit between snow albedo radiative forcing change and temperature increase during the study period, Chen et al. calculated snow albedo feedback over the North Hemisphere to be 0.17 ± 0.008 W∙m⁻²∙K⁻¹ [57]. By multiplying the sensitivity of end date of snow cover to temperature by difference in black surface albedo before and after snowmelt, the result is 0.19 ± 0.17 W∙m⁻²∙K⁻¹, according to Peng et al. [58]. With the method of block bootstrap test, snow albedo feedback of this study is 0.18 ± 0.08 W∙m⁻²∙°C⁻¹. Despite the different methods and datasets applied in each work, results of the three studies agree well.

Snow albedo radiative forcing is also compared, as it’s a key parameter of feedback and offers valuable spatial information. Both Flanner and coworkers’ work and our work show similar snow albedo radiative forcing pattern in terms of spatial distribution and temporal variation [33]. Large snow albedo radiative forcing is found in the northern part of Eurasia, the mid-high latitude of North America and the Tibetan Plateau in both works, though the Tibetan Plateau exhibits even larger values in our work. In addition, seasonal cycle of snow albedo radiative forcing in both works show a single peak (peaks in April), with the largest value in spring months. Chen and coworkers’ work, however, shows a relatively weaker consistency with us [57]: smaller value in the Tibetan Plateau as compared with Flanner et al. and our work, and a clear tendency of snow albedo radiative forcing with latitude, i.e., the higher latitudes exhibit larger radiative forcing.

While different surface albedo kernels used in these studies have proven to vary only a little [33,56], differences among the results of snow albedo radiative forcing are considered mainly due to albedo change caused by snow cover change ($\frac{\partial \alpha }{\partial S}(t,r)$). In Chen and coworkers’ work, surface albedo instead of snow albedo data was used, thus the influence of vegetation change was imported to snow albedo change [57]. In addition, according to the conclusions of Singh et al., the coarse resolution of snow data (1° and monthly [33]) is likely to result in an overestimate of snow albedo radiative forcing [53].

4.2. Comparison with Model-Based Studies

Based on the agreement among observation-based studies, snow albedo feedbacks of our work are compared with those from the 25 models that participated in CMIP5 (Qu and Hall [56], Table 1 of their paper).

Snow albedo feedback over the NEL estimated from the 25 models ranges from 0.18 to 0.78 W∙m⁻²∙K⁻¹, with the ensemble mean of 0.42 ± 0.15 W∙m⁻²∙K⁻¹. Our result is 0.18 ± 0.08 W∙m⁻²∙K⁻¹. The global snow albedo feedback of the 25 models ranges from 0.03 to 0.16 W∙m⁻²∙K⁻¹, with ensemble mean of 0.08 ± 0.03 W∙m⁻²∙K⁻¹, and, in our case, it is 0.04 ± 0.02 W∙m⁻²∙K⁻¹. In general, results of our work fall near the lower bound of the 25 models. The ensemble means of both hemispheric and global snow albedo feedbacks estimated by Qu and Hall [56] are larger than our results and the other partially observation-based results mentioned above. This might indicate an overall overestimation of snow albedo feedback by most of the 25 models. For the purpose of offering detailed information for model optimism, the possible source of discrepancies between this study and the 25 models is discussed below.

According to Qu and Hall, snow albedo feedback is mainly determined by two terms: one represents the variations in planetary albedo with surface albedo ($\frac{\partial {\alpha }_P}{\partial {\alpha }_S}$), and the other is the change in surface albedo associated with a 1 °C increase in T_s ($\frac{∆{\alpha }_S}{∆T_S}$) [56]. The spread of $\frac{\partial {\alpha }_P}{\partial {\alpha }_S}$is mainly determined by the surface albedo kernel, because this term is calculated as the ratio of surface albedo kernel to the TOA shortwave radiation. As the kernels used in models are also used in this study, we consider this coefficient is of little contribution to the difference. As a result, the major source of the discrepancy would be originated from $\frac{∆{\alpha }_S}{∆T_S}$. Specifically, increase in the area-averaged surface air temperature ($∆Ts$) is relatively consistent both among the models and between model simulations and observations, thus surface albedo change ($∆{\alpha }_S$) could be the largest source of discrepancy between our work and the 25 models. Surface albedo change, which is mainly determined by albedo contrast between snow-covered and snow-free surface, is also considered to be most uncertain (threefold spread in CMIP3 and persists to be large in CMIP5) among models. However, as the spread of surface albedo change between snow-covered and snow-free surface is not directly analyzed in Qu and Hall’s work [56], further investigations are still required.

Despite the fact that there are other factors contributing to the discrepancies, surface albedo change is considered as the predominate difference and possibly being overestimated by most of the 25 models. Therefore, model parameterization should specifically focus on this factor.

Moreover, according to Hall and Qu, the surface albedo decrease associated with loss of snow cover, rather than the reduction in snow albedo due to snow metamorphosis is more important in the determination of snow albedo feedback [41]. Thus, in turn, the significance of albedo constraint by snow cover data of high spatial and temporal resolutions in this study is strengthened. The information of snow cover change, as well as the constrained albedo data, can be guidance for the model parameterizations.

4.3. Strengths and Limitations

The strengths of this study compared with the previous can be summarized as the following:

• Instead of a combination assessment of snow and ice albedo feedback, snow albedo feedback is examined exclusively in this study, thus the contribution of snow albedo feedback and its uncertainty to the surface albedo feedback can be independently achieved.
• Satellite-based MODIS snow cover data of high spatial resolution (0.05°) is used to constrain and determine the areas of snow albedo radiative forcing and feedback. In this study, when snow cover data were included in the calculation, the snow albedo radiative forcing decreased as much as 27.63% (not shown), compared to when only albedo contrast data and surface albedo kernel data were used for calculation. Meanwhile, it offers accurate information on constraining surface albedo decrease associated with loss of snow cover of models on spatial aspect.
• High temporal resolution (daily) of snow cover and albedo data offers detailed information of snow cover change, which is relevant because changing snow cover occurs rapidly. In particular, the snow melt in mid-latitudes generally lasts for less than one month, especially in spring. For example, the Tibetan Plateau, one of the most intensive and important snow albedo feedback areas, has the largest inconsistency in snow albedo radiative forcing according to the comparisons previously (Section 4.1), since rapid snow accumulation and ablation processes in spring always last for less than one week. Thus, the monthly mean data, which are commonly used in other studies [35,43,59], would easily smooth the snow and feedback processes. In addition, surface albedo decrease associated with loss of snow cover is estimated precisely on temporal scale, which offers guidance for model optimization.
• The block bootstrap test is used effectively to reduce the uncertainty of snow albedo feedback. Considering the fact that most observation-based studies are short in time duration, a simple linear regression [28,35,43] to compute snow albedo feedback and its confidence limits would probably give misleading results due to the random variations of variables. By enlarging the sample amount and enhancing the number of tests (10,000 times in our case), the block bootstrap test should have obtained more reliable results.

Our study has some limitations that could lead to the uncertainties of the results:

• Compared to model simulations, the available observational data of only 14 years would be potentially biased by internal climate variation.
• Data from multiple sources with different spatial and temporal resolutions are applied to our work. Therefore, the processes of unifying their spatial and temporal resolutions, i.e., interpolation and resampling, would add errors to the spatial distribution of the results. Meanwhile, different temporal resolutions between daily datasets (fractional snow cover and albedo contrast) and monthly dataset (the radiative kernel) would also add uncertainty in temporal variation of the results.
• The import of snow cover data is intended to constrain albedo data, quantifying a more precise area of snow albedo radiative forcing and feedback. However, due to different data sources and different computation methods, discrepancies are imported as well, and would probably be the major source of uncertainty.
• Long-term observational data are not available for feedback study at present, and, at the same time, intermodel spread cannot be constrained effectively. Therefore, the best effort we can make is probably observing short-term variations and comparing the results with those from climate models [22,60]. Even though there has not been substantial progress in using observation results to constrain model simulations directly, we believe that, by using improved observational datasets and methods, observation-based results would help in better understanding the origin of intermodel differences, as well as the assessment of reliability of different model simulations. Finally, the goal is to get better description of feedback processes and finer estimation of feedbacks, more accurate ECS, and better projections of future climate [2,16].

5. Conclusions

We quantified snow albedo feedback and calculated snow albedo radiative forcing during the period of 2003–2016, with remote sensing, atmospheric reanalysis and radiative kernel data. The results were then compared with those from partially observation-based calculations and model-based estimations.

The results suggest that, excluding Antarctica, both snow cover and snow albedo radiative forcing are the largest in landmasses north of 30° N. Snow albedo radiative forcing peaks in spring, due to strong insolation and large snow cover extent. Snow albedo feedback over the NEL is estimated to be 0.18 ± 0.08 W∙m⁻²∙°C⁻¹ and the global mean is 0.04 ± 0.02 W∙m⁻²∙°C⁻¹. Results were compared with partially observation-based studies first. The regional mean snow albedo feedbacks were consistent, while the spatial pattern of snow albedo radiative forcing was different. The differences probably originated from snow albedo data and albedo data, as well as resolutions of data. Compared to the 25 models that participated in CMIP5, a general overestimation of models was found, mainly due to the overestimation of surface albedo variation between snow-covered and snow-free surface in the models. Surface albedo change is also with the largest spread among snow albedo feedback determinants. Therefore, model parameterization should specifically focus on the constraint of this factor. Meanwhile, remotely sensed snow cover data with high spatial and temporal resolutions, the constrained albedo contrast data between snow-covered and snow-free surface, as well as the resulted high resolution of snow albedo radiative forcing in this study, can offer valuable information for model parameterization.