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].

**Figure 1.** Spatial distribution of monthly surface albedo radiative kernel under all sky (2.8◦ on average resampled to 0.25◦ cell size) (these data are based on NCAR's monthly surface albedo radiative kernel [47]).
