*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–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\_{\mathbf{s}}(t,\mathsf{R}) \; , \; \begin{cases} \begin{array}{c} 0, \\ \frac{\int\_{\mathbb{R}} \frac{\partial \mathbf{u}}{\partial \mathbf{S}}(t,r) \times \frac{\partial H}{\partial \mathbf{n}}(t,r) \times dA(r)}{A(\mathsf{R})} \end{array} , \; S(t,r) > 0 \end{cases} \tag{1}$$

Here, *Gs(t,R)* is snow albedo radiative forcing that describes the instantaneous influence of snow cover on TOA energy budget, with unit of *<sup>W</sup>*·*m*−*<sup>2</sup>* [33]. *Gs(t,R)* is a function of time *<sup>t</sup>* and study area *<sup>R</sup>*, 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., *Gs(t,R)* = 0. Only pixels with a fractional snow cover larger than 0 contribute to the snow albedo radiative forcing. The magnitude of *Gs(t,R)* depends on the albedo contrast between snow-covered and snow-free circumstances, as well as the radiative kernel data. Specifically, *∂α <sup>∂</sup><sup>S</sup>* (*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]). *<sup>∂</sup><sup>H</sup> ∂α* (*t*,*r*) is the radiative change against the coincident change of albedo, namely the surface albedo radiative kernel. Note that *∂α <sup>∂</sup><sup>S</sup>* (*t*,*r*) and *S(t,r)*, *<sup>∂</sup><sup>H</sup> ∂α* (*t*,*r*) and *α(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]:

$$\lambda = \frac{\text{Gs}(t, R)}{T(t, R)}\tag{2}$$

where <sup>λ</sup> is the feedback parameter, here considered as snow albedo feedback, with unit of <sup>W</sup>· <sup>m</sup>−2· ◦C<sup>−</sup>1. Δ*Gs*(*t*,*R*) and Δ*T*(*t*,*R*) are monthly anomalies of *Gs*(*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, *Gs*(*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 *Gs*(*t*,*R*) and *T*(*t*,*R*) were area weighted [34] during the calculation process.

The coefficient of least square fit of Δ*Gs*(*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 Δ*Gs*(*t*,*R*) and Δ*T*(*t*,*R*), the result may be subjected to substantial uncertainty due to random variations of *Gs*(*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 Δ*Gs*(*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.

**Figure 2.** Spatial distribution of mean fractional snow cover during 2003–2016.

**Figure 3.** Fractional snow cover over the North Hemisphere Extratropical Land during 2003–2016: (**a**) monthly and annual mean; and (**b**) climatological monthly mean.
