Contribution of Changes in Snow Cover Extent to Shortwave Radiation Perturbations at the Top of the Atmosphere over the Northern Hemisphere during 2000–2019

: Snow-induced radiative forcing (S n RF), deﬁned as the instantaneous perturbation of the Earth’s shortwave radiation at the top of the atmosphere (TOA), results from variations in the terrestrial snow cover extent (SCE), and is critical for the regulation of the Earth’s energy budget. However, with the growing seasonal divergence of SCE over the Northern Hemisphere (NH) in the past two decades, novel insights pertaining to S n RF are lacking. Consequently, the contribution of S n RF to TOA shortwave radiation anomalies still remains unclear. Utilizing the latest datasets of snow cover, surface albedo, and albedo radiative kernels, this study investigated the distribution of SnRF over the NH and explored its changes from 2000 to 2019. The 20-year averaged annual mean S n RF in the NH was − 1.13 ± 0.05 W m − 2 , with a weakening trend of 0.0047 Wm − 2 yr − 1 ( p < 0.01) during 2000–2019, indicating that an extra 0.094 W m − 2 of shortwave radiation was absorbed by the Earth climate system. Moreover, changes in S n RF were highly correlated with satellite-observed TOA shortwave ﬂux anomalies ( r = 0.79, p < 0.05) during 2000–2019. Additionally, a detailed contribution analysis revealed that the S n RF in snow accumulation months, from March to May, accounted for 58.10% of the annual mean S n RF variability across the NH. These results can assist in providing a better understanding of the role of snow cover in Earth’s climate system in the context of climate change. Although the rapid SCE decline over the NH has a hiatus for the period during 2000–2019, S n RF continues to follow a weakening trend. Therefore, this should be taken into consideration in current climate change models and future climate projections.


Introduction
Snow cover is an integral component of the cryosphere and represents one of the three essential climate variables (ECVs) for snow in the Global Climate Observing System (GCOS) [1], which plays a crucial role in the Earth's climate system by regulating the surface energy budget [2][3][4]. As a result of the high reflectivity of snow cover, the presence or absence of snow cover plays a key role in the regulation of the heating and cooling patterns of the Earth's surface, more than any other land surface category. For example, as the surface air temperature increases, the snow cover extent (SCE) retreats and the surface albedo decreases. Consequently, the land surface reflects less solar radiation back to space and retains more energy. The reduced reflected shortwave radiation at the top of the atmosphere (TOA) is known as snow-induced radiative forcing (S n RF) [2,5,6], which is critical for regulating Earth's energy budget.
As reported by Flanner et al. [2] and Chen et al. [6], S n RF (W m −2 ) can be calculated by combining observations of SCE, surface albedo contrast (∆α) between snow-covered and snow-free surfaces, and TOA flux variations with surface albedo change (∂F/∂α). The time (t) dependence of S n RF, within a region R (here, the Northern Hemisphere; NH) of area A, can be expressed as Equation (1) [2,5]: where S (t, r) is the range of snow cover over the NH; α (t, r) is the surface albedo; ∂α/∂S(t, r) is the rate of variation of the surface albedo with changes in snow cover; and ∂F/∂α(t, r) is the response of the TOA net shortwave radiation anomalies to surface albedo changes. According to Flanner et al. [2] and Chen et al. [6], we assumed that ∂α/∂S and ∂F/∂α were constant for the temporally and spatially varying snow cover S and surface albedo α, respectively. Furthermore, ∂α/∂S could be replaced by the surface albedo contrast ∆α driven by the snow cover anomaly, and ∂F/∂α could be substituted by albedo radiative kernels [2,5,6], wherein the albedo radiative kernels are expressed as the TOA shortwave flux anomalies associated with a 1% change in surface albedo. The contribution of SCE to the Earth's energy budget has received substantial attention in recent decades. Several studies have estimated S n RF at regional and global scales. For example, S n RF in snow accumulation and melting seasons over the NH was observed to be 0.01 and 0.12 W m −2 during 1982-2013, respectively [3]. However, SCE reduction was observed to result in the absorption of approximately 0.22 W m −2 extra shortwave radiation in the NH during 1979-2008, which is nearly equal to the radiative forcing caused by Arctic sea ice shrinkage [2]. Similarly, driven by advanced snow end dates, S n RF in the NH was weakened by approximately 0.16 W m −2 from 2001 to 2013 [5]. As per the representative concentration pathway (RCP) 8.5 scenario, a reduction in SCE over the NH may lead to an extra absorption of 0.47-0.60 W m −2 of shortwave radiation in the 21st century [7]. The above-mentioned studies have proved the contribution of SCE to Earth's radiation budget. However, owing to the uncertainty of input datasets, with a poor quality of the surface albedo dataset, incomplete spatial coverage of snow cover dataset, and short time span of direct TOA radiation observations, significant differences exist in the S n RF estimations from climate models, reanalysis datasets, satellite products, and ground observations. For example, estimates of the S n RF over the NH from the moderate resolution imaging spectroradiometer (MODIS) were approximately 18% lower than the comparable estimates from the coarse-resolution advanced very high resolution radiometer (AVHRR) during 2001-2008 [8]. Moreover, S n RF is a key parameter for the calculation of snowalbedo feedback (SAF), which may have significant effects on climate projection studies. For example, SAF exhibits a five-fold spread across the coupled model intercomparison project phase 5 (CMIP5) models, ranging from 0.03 to 0.16 W m −2 K −1 [9]. In addition, SAF from ground observations is apparently higher than the comparable results from the reanalysis dataset [10]. Therefore, an accurate estimation of S n RF is required for measuring the influence of snow cover on Earth's radiation budget [2,8], reducing the bias of SAF calculation, and improving climate projections [7,9]. S n RF is determined by several factors, including SCE change, surface albedo anomalies induced by snow cover anomalies (which rely on snow-covered and snow-free surface albedo), solar insolation, and atmospheric states, which determine the propagation of surface albedo changes to the TOA radiation flux [2]. Compared with other variables, snow cover change is the most important driving factor of S n RF. Without snow cover, the ground absorbs approximately four to six times more energy from the sun [11,12]. Instead of a continuous rapid decline in SCE, a recent study reported a notable seasonal divergence of terrestrial SCE across the NH over recent decades [13]. Considering the importance of snow cover for estimating S n RF, quantifying and understanding its spatial and temporal variations is essential for determining its contribution to the Earth's energy budget.
Optimizing the snow cover and surface albedo inputs utilizing satellite observations is a key approach for improving the accuracy of S n RF [14]. To quantify S n RF and explore its contribution to TOA shortwave radiation anomalies driven by the latest SCE features in the NH, a novel estimation of S n RF using satellite observations at the best achievable spatial and temporal resolutions is urgently required. Therefore, the objective of the present study was two-fold: (1) to quantify the changes in SCE and resulting S n RF; (2) to determine the contribution of S n RF to TOA energy budget anomalies. To achieve these objectives, we first investigated the latest SCE anomalies over the NH. Furthermore, we chose the surface albedo dataset with the best performance resulting from a cross-comparison of SCE changes. Finally, we calculated the S n RF based on the selected surface albedo and explored its spatial distribution and spatial-temporal changes, as well as its influence on the Earth's energy budget.

Materials
To meet the requirements of the S n RF calculation, three variables-snow cover extent, surface albedo, and albedo radiative kernels-were employed in this study.

Snow Cover Extent Dataset
To detect SCE anomalies in the NH and to implement them as input parameters for the estimation of S n RF, a monthly snow chart dataset, calculated from the NH SCE CDR v01r01 [15] and the snow cover range from the MODIS/Terra monthly SCE in 0.05 • climate modeling grid (CMG) cells (MOD10CM) [16], were employed in this study.

NHSCE Snow Charts
Monthly snow charts were obtained from the NH SCE CDR v01r01 and utilized to determine the long-term SCE anomaly in the NH from 1972 to 2020. These snow chart data are widely implemented in large-scale SCE anomaly studies, owing to their consistency and long time span [15].

MOD10CM
MOD10CM is an aggregation of MODIS/Terra Snow Cover Daily L3 Global 0.05 • CMG (MOD10C1) products [17], providing the monthly mean snow cover percentage in each 0.05 • grid cell from February 2000 to the present. The monthly composite of daily MOD10C1 is essential in climate change studies because of the persistent cloudiness and limited spatial coverage of MODIS daily snow cover products, particularly at high latitudes [18]. In this study, MOD10CM was utilized as the snow-cover range in the S n RF calculation.
Compared with other snow cover products, MOD10CM is the only consistent, objective snow cover estimate retrieved from optical satellite observations in the past 20 years with a relatively fine spatial resolution. The 20-year averaged snow cover percentage in the NH calculated from MOD10CM during 2000-2019 is displayed in Figure 1. To match the limited spatial coverage of the surface albedo dataset, high latitudes above 75 • N were excluded in this study.

Land Surface Albedo Dataset
The surface albedo, defined as the ratio of the reflected radiation flux to the incoming radiation flux, is one of the most important parameters for calculating S n RF, as it indicates the amount of solar radiation reflected back to space because of snow cover [2,5,7,8,19], and has been further designed as one of the GCOS ECVs [20].

Land Surface Albedo Dataset
The surface albedo, defined as the ratio of the reflected radiation flux to the incoming radiation flux, is one of the most important parameters for calculating SnRF, as it indicates the amount of solar radiation reflected back to space because of snow cover [2,5,7,8,19], and has been further designed as one of the GCOS ECVs [20].
To quantify ∆α, required for estimating SnRF, the broadband directional surface albedo (ALBB-DH) produced by the Copernicus Climate Change Service (C3S) at a spatial resolution of 1 km during 2000-2019 [21] and the MODIS Terra and Aqua combined bidirectional reflectance distribution function (MCD43GF) at a spatial resolution of 1 km during 2001-2017 [22] were implemented in this study. C3S ALBH-DH is the land surface albedo, and MCD43GF is the snow-free surface albedo. Thus, we obtained yearly and monthly ∆α by subtracting MCD43GF from C3S ALBH-DH in the corresponding period.
The broadband directional surface albedo (ALBB-DH) produced by the Copernicus Climate Change Service (C3S) [21] is the best option for surface albedo inputs in SnRF estimation because of its fine spatial resolution (1 km) and complete spatial coverage. However, the C3S ALBH-DH is a newly released dataset. To verify the fitness of C3S ALBH-DH for SnRF estimation, we employed two widely utilized surface albedo datasets, including the global land surface satellite (GLASS) [23] and the Cloud, Albedo, and Surface Radiation datasets from AVHRR data Edition 2 (CLARA-SAL-A2) [24,25] for crosscomparison analysis.

C3S ALBH-DH Surface Albedo Dataset
The C3S ALBH-DH is a global 10-daily gridded surface albedo dataset with a spatial resolution of 1 km covering 1998-2020, calculated from the Systeme Probatoire d'Observation de la Terre-VEGETATION (SPOT-VGT) and the Project for On-Board Autonomy-Vegetation (PROBA-V) [21]. In this study, the integrated black-sky surface albedo computed over the total spectrum (0.4-4.0 µm) from 2000 to 2019 was utilized for analysis. Compared with other surface albedo datasets, such as GLASS, CLARA-SAL-A2, MODIS, the Advanced Polar Pathfinder (APP-X), and the International Satellite Cloud Climatology Project (ISCCP), C3S ALBH-DH is ideal for continental-scale SnRF estimation because of its fine spatial resolution, long time span, and spatial completeness.

GLASS Land Surface Albedo Dataset
The GLASS surface albedo dataset was developed by integrating AVHRR and MODIS version 6 radiance data through a temporal filter scheme [23]. The GLASS surface albedo dataset provides surface albedo at total shortwave, visible, and near-IR spectral To quantify ∆α, required for estimating S n RF, the broadband directional surface albedo (ALBB-DH) produced by the Copernicus Climate Change Service (C3S) at a spatial resolution of 1 km during 2000-2019 [21] and the MODIS Terra and Aqua combined bidirectional reflectance distribution function (MCD43GF) at a spatial resolution of 1 km during 2001-2017 [22] were implemented in this study. C3S ALBH-DH is the land surface albedo, and MCD43GF is the snow-free surface albedo. Thus, we obtained yearly and monthly ∆α by subtracting MCD43GF from C3S ALBH-DH in the corresponding period.
The broadband directional surface albedo (ALBB-DH) produced by the Copernicus Climate Change Service (C3S) [21] is the best option for surface albedo inputs in S n RF estimation because of its fine spatial resolution (1 km) and complete spatial coverage. However, the C3S ALBH-DH is a newly released dataset. To verify the fitness of C3S ALBH-DH for S n RF estimation, we employed two widely utilized surface albedo datasets, including the global land surface satellite (GLASS) [23] and the Cloud, Albedo, and Surface Radiation datasets from AVHRR data Edition 2 (CLARA-SAL-A2) [24,25] for cross-comparison analysis.

C3S ALBH-DH Surface Albedo Dataset
The C3S ALBH-DH is a global 10-daily gridded surface albedo dataset with a spatial resolution of 1 km covering 1998-2020, calculated from the Systeme Probatoire d'Observation de la Terre-VEGETATION (SPOT-VGT) and the Project for On-Board Autonomy-Vegetation (PROBA-V) [21]. In this study, the integrated black-sky surface albedo computed over the total spectrum (0.4-4.0 µm) from 2000 to 2019 was utilized for analysis. Compared with other surface albedo datasets, such as GLASS, CLARA-SAL-A2, MODIS, the Advanced Polar Pathfinder (APP-X), and the International Satellite Cloud Climatology Project (ISCCP), C3S ALBH-DH is ideal for continental-scale S n RF estimation because of its fine spatial resolution, long time span, and spatial completeness.

GLASS Land Surface Albedo Dataset
The GLASS surface albedo dataset was developed by integrating AVHRR and MODIS version 6 radiance data through a temporal filter scheme [23]. The GLASS surface albedo dataset provides surface albedo at total shortwave, visible, and near-IR spectral ranges under actual atmospheric conditions. To meet the objectives of the present study, 8-day GLASS total shortwave black-sky surface albedo data at a spatial resolution of 0.05 • , from 2000 to 2019, were utilized in this study. The GLASS surface albedo dataset has previously been used to quantify the albedo-induced radiative forcing of snow melting over Greenland [26] and snow cover phenology changes over the NH [3,5]. The CLARA-SAL-A2 surface albedo dataset was generated based on a homogenized AVHRR radiance time series and was created using algorithms to derive surface albedo for different land use types, including snow, sea ice, open water, and vegetation. Currently, CLARA-SAL-A2 is the only available long-time-span surface albedo product that is completely dependent on AVHRR imagery [24,25]. For the purpose of the present study, the monthly means of CLARA-SAL-A2 at a spatial resolution of 0.25 • from 2000 to 2019 were employed in the analysis.

MCD43GF Snow-Free Surface Albedo Dataset
To obtain the climatology of monthly mean snow-free surface albedo required in S n RF estimation, a global gap-filled snow-free surface albedo dataset retrieved from the MODIS Terra and Aqua combined bidirectional reflectance distribution function (MCD43GF) at a spatial resolution of 1 km during 2001-2017 [22] was utilized in this study. MCD43GF provides spatially complete surface albedo values after removing the presence of snow cover, which are computed by MODIS spectral bands 1-7, as well as shortwave infrared, visible, and near-infrared bands. For the present study, the daily broadband snow-free black-sky surface albedo data from 2001 to 2017 were used.

Albedo Radiative Kernels Datasets
The albedo radiative kernel, expressed as the TOA shortwave flux anomaly associated with a 1% change in surface albedo, is a key input variable for estimating S n RF. Compared with analytical models, the radiative kernel approach has significant advantages in separating the radiative responses to land surface albedo induced by snow cover changes from other climate parameters, enhancing the accuracy for capturing the functional dependence of planetary albedo on surface albedo [9].
To reduce uncertainties in S n RF estimation, three albedo radiative kernels, the National Center For Atmospheric Research (NCAR) Community Atmosphere Model 3 (CAM3) gridded at a horizontal resolution of 2.81 • [27], the Geophysical Fluid Dynamics Laboratory (GFDL) Atmospheric Model 2 (AM2) gridded at a horizontal resolution of 2.50 • [28], and the sixth generation atmospheric general circulation model (ECHAM6) gridded at a horizontal resolution of 1.875 • [29], were used in this study. The intra-annual variability of the three albedo radiative kernels is displayed in Figure 2.

Data Preparation
A summary of the datasets implemented in this study is presented in Table 1. To cater to the requirements of the SnRF calculation, the snow cover range, land surface albedo, and albedo radiative kernels were regridded at a spatial resolution of 0.05° using the gdalwarp reprojection and warping utility (https://gdal.org/programs/gdalwarp.html, Figure 2. Three albedo radiative kernels used in this study.

TOA Radiation Budget Dataset
To explore the contribution of S n RF to TOA energy budget anomalies, the clear-sky upward TOA shortwave (0.2-5 µm) flux, retrieved from the Clouds and Earth's Radiant Energy Systems (CERES) Energy Balanced and Filled (EBAF) products [30] at a spatial resolution of 1.0 • , from 2000 to 2019, was implemented. The CERES EBAF clear-sky shortwave flux is spatially complete, and is inferred from both CERES and MODIS measurements conducted monthly [31].

Data Preparation
A summary of the datasets implemented in this study is presented in Table 1. To cater to the requirements of the S n RF calculation, the snow cover range, land surface albedo, and albedo radiative kernels were regridded at a spatial resolution of 0.05 • using the gdalwarp reprojection and warping utility (https://gdal.org/programs/gdalwarp.html, acessed on 28 November 2021). For datasets with a spatial resolution coarser than 0.05 • , we used "cubicspline" in the resampling process; for datasets with spatial resolution finer than 0.05 • , we utilized "average" in the resampling process, which computes the weighted average of all contributing pixels (excluding NODATA) at 0.05 • . Moreover, to match the temporal resolution (monthly) of albedo radiative kernels, both the 10-day C3S-ALBH-DH and daily MCD43GF were reproduced to generate monthly mean values.

Z-Score Normalization
Z-score normalization was implemented to normalize the different variables in the crosscomparison analysis. For variable X, the Z-score (Z x ) was calculated using Equation (2) [32].
where µ x is the average value of X and δ x is the standard deviation of X.

Linear Slope
The linear slope measures the rate of change in the dependent variable as the independent variable changes, which can be calculated using Equation (3): where ∆y is the change in dependent variable y and ∆x is the change in independent variable x.

Pearson Correlation Coefficient
The Pearson correlation coefficient (r) was utilized to quantify the degree of correlation between two variables. For variables X and Y, r was calculated using Equation (4): where cov(x, y) is the covariance of variable x and variable y; δ x is the standard deviation of variable x; and δ y is the standard deviation of variable y.

Relative Contribution Calculation
We assumed that the variability in the annual mean S n RF was attributed to S n RF anomalies in each month. Furthermore, regression analysis was utilized to quantify the contribution of S n RF in each month to the annual mean S n RF variability during 2000-2019 using Equation (5): where S n RF z is the Z-score of the normalized annual mean S n RF during 2000-2019; S n RF iz is the Z-score normalized S n RF in month i (i = 1, 12); β i is the regression coefficient for S n RF iz ; and ε is the residual error. Furthermore, the relative contribution (C i ) of S n RF iz to S n RF z was confirmed using Equation (6): where S n RF iz_slope is the linear slope of S n RF iz , which can be calculated using Equation (3). The regression analysis is an effective approach for solving collinearity problems in contribution analysis [33] and has been applied in gross primary production trend analysis in the Three-North region of China [34], as well as for determining the interannual evapotranspiration sensitivity to climate change across the Yellow River basin [35].

Results
To explore the contribution of changes in SCE to TOA shortwave radiation perturbation in the NH, we first investigated the SCE anomalies during 2000-2019. Further, we verified the performance of the C3S surface albedo through a cross-comparison with the snow cover dataset. Finally, we calculated the variabilities of the annual mean S n RF and identified the dominant month that controlled the annual mean S n RF.

Characteristic of SCE Variability in the NH during 2000-2019
The Z-score anomalies of annual mean SCE in the NH, including North America (NA) and Eurasia (EU), are displayed in Figure 3a. The climatology of SCE in each month and the associated changes during 2000-2019 are exhibited in Figure 3b,c, respectively. Changes were calculated by the monthly mean SCE during 2000-2019 minus the comparable values in the referencing period, 1972-2000.
As displayed in Figure 3a, the well-documented rapid SCE decline over the NH since the late 1960s [36][37][38] has experienced a hiatus in recent years, especially during 2000-2019, both in NA and EU. The SCE in the NH declined at a speed of −0.52 × 10 6 km 2 per decade (p < 0.05) during 1972-2000, but displayed negligible changes during 2000-2019. The ignorable SCE anomalies from 2000 to 2019 disagree with climate projections, wherein the SCE will continue to shrink in the future [37,39,40]. rable values in the referencing period, 1972-2000.
As displayed in Figure 3a, the well-documented rapid SCE decline over the NH since the late 1960s [36][37][38] has experienced a hiatus in recent years, especially during 2000-2019, both in NA and EU. The SCE in the NH declined at a speed of -0.52 × 10 6 km 2 per decade (p < 0.05) during 1972-2000, but displayed negligible changes during 2000-2019. The ignorable SCE anomalies from 2000 to 2019 disagree with climate projections, wherein the SCE will continue to shrink in the future [37,39,40]. The seasonal variability of SCE in the NH is presented in Figure 3b. The maximum SCE in the NH occurred in January (47.78 ± 1.59 10 6 km 2 ) during 2000-2019. Meanwhile, The seasonal variability of SCE in the NH is presented in Figure 3b. The maximum SCE in the NH occurred in January (47.78 ± 1.59 10 6 km 2 ) during 2000-2019. Meanwhile, the minimum SCE in the NH occurred in August (2.59 ± 0.63 10 6 km 2 ) during the same period. Based on the seasonal cycle of SCE, and published studies, such as those by Chen et al. [13] and Brutel-Vuilmet et al. [41], the months from March to August were defined as the snow melting season and the months from September to February were the snow accumulation season. Detailed SCE changes per month further revealed that, in comparison with 1972-2000, SCE anomalies between accumulation and melting seasons were markedly different during 2000-2019 (Figure 3c). In the melting season, SCE in the NH declined from 1.86 × 10 6 km 2 in June to 0.19 × 10 6 km 2 in March. Conversely, SCE increased in the accumulation season over the NH during the same period, ranging from 0.26 × 10 6 km 2 in September to 2.35 × 10 6 km 2 in October. The distinct SCE anomalies in each month indicate that, instead of a rapid SCE decline, the intra-annual SCE variability is a new characteristic of snow cover in recent years.  (Figure 3c). In the melting season, SCE in the NH declined from 1.86 × 10 6 km 2 in June to 0.19 × 10 6 km 2 in March. Conversely, SCE increased in the accumulation season over the NH during the same period, ranging from 0.26 × 10 6 km 2 in September to 2.35 × 10 6 km 2 in October. The distinct SCE anomalies in each month indicate that, instead of a rapid SCE decline, the intra-annual SCE variability is a new characteristic of snow cover in recent years.

Validation of C3S Surface Albedo Dataset in SnRF Calculation
A cross-comparison was conducted for SCE and surface albedo anomalies between March and June, which presented the minimum and maximum negative changes in SCE over the NH during 2000-2019.  The Z-scores of SCE for March and surface albedo from C3S ALBH-DH, GLASS, and CLARA-SAL-A2, and linear scatter plots between the Z-scores of SCE for March and surface albedo during 2000-2019 are displayed in Figure 4e,f, respectively. The performance of surface albedo is largely determined by its spatial resolution. Compared with GLASS (r = 0.72, p < 0.05) and CLARA-SAL-A2 (r = 0.64, p < 0.05), C3S ALBH-DH (r = 0.76, p < 0.05) displayed an improved linear correlation coefficient with SCE during 2000-2019 (Figure 4f). This indicates that C3S ALBH-DH skillfully captured the SCE variability and performed better than GLASS and CLARA-SAL-A2 during 2000-2019. Moreover, uncertainty in snow cover data [36,37] and changes in snow morphology [44] may reduce its correlation with surface albedo. However, C3S ALBH-DH still captured SCE anomalies in March during 2000-2019.

Performance of C3S-ALBH-DH in June during 2000-2019
The 20-year average snow cover percentage in June and the surface albedo from C3S ALBH-DH, CLARA-SAL-A2, and GLASS from 2000 to 2019 are exhibited in Figure 5.

Estimation of SnRF over the NH
Applying the albedo radiative kernel approach using Equation (1)  The SCE in June decreased nearly twice as rapidly as the widely acknowledged sea ice extent shrinking in September during 1979-2011 [38] and exhibited a sustained rapid decrease during 2000-2019 (Figure 3c). To verify the ability of C3S ALBH-DH at capturing the rapid SCE decline in June, we first investigated the spatial distribution of the snow cover percentage and surface albedo in June. Further, we calculated the correlation coefficient between changes in the SCE and surface albedo in June during 2000-2019.
As shown in Figure 5a, snow cover in June was mostly distributed at high latitudes in NA and EU. Meanwhile, the surface albedo from C3S ALBH-DH (Figure 5b)

Estimation of S n RF over the NH
Applying the albedo radiative kernel approach using Equation (1), we first calculated the S n RF during 2000-2019. Further, we compared the S n RF estimates with direct TOA shortwave flux observations and explored the relative contribution from each month to the annual mean S n RF from 2000 to 2019.

Estimated S n RF over the NH during 2000-2019
The spatial distributions of the 20-year annual mean S n RF in the NH under clear-sky conditions, and the changes that occurred during 2000-2019, are displayed in Figure 6a  The 20-year annual mean S n RF in NH was estimated to be −1.13 (±0.05) W m −2 from 2000 to 2019 (Figure 6a,c), with high values of S n RF distributed at high latitudes in NA and EU, the Pamir Mountains, and East Siberia. Moreover, the change in S n RF over the NH was estimated to be −0.1033 (± 0.0081) W m −2 from 2000 to 2019. The absolute value of S n RF is decreasing, indicating a generally receding trend of radiative cooling in the NH from 2000 to 2019. Meanwhile, the minima of annual mean S n RF occurred in 2009 and 2011, coinciding with the large-scale El Niño and Arctic Oscillation events, particularly in 2009 [42]. Conversely, the peaks occurred in 2015 and 2019, corresponding to La Niña and extremely high temperatures [45].
The intra-annual cycles of S n RF in the NH during 2000-2019 indicated that the S n RF peaks at approximately −3.83 (± 0.14) W m −2 in April, when both solar irradiance and SCE over the NH are high. Compared with the S n RF in the NH during 1979-2008 [2], the peak of S n RF was much lower. The detailed S n RF anomalies in each month revealed that S n RF significantly decreased from February to June, indicating a weakened radiative cooling effect over the NH during these months in 2000-2019. The most significant S n RF decline occurred in May at 0.3681 W m −2 (p < 0.01) from 2000 to 2019. Meanwhile, S n RF weakened by 0.1714 W m −2 (p < 0.01) in June and 0.1559 W m −2 (p < 0.05) in March. Conversely, S n RF exhibits negligible changes in snow accumulation months, especially in October, November, January, and February. The contrasting S n RF changes in snow accumulation and melting seasons may lead to unbalanced intra-annual shortwave radiative distribution, which should be considered in climate change projections.

Direct Estimates of TOA Shortwave Flux Anomalies
The 20-year averaged annual mean TOA shortwave flux and changes calculated from CERES EBAF observations from 2000 to 2019 are displayed in Figure 7a,b. The variability of annual mean TOA upward shortwave flux and its linear correlation with S n RF are exhibited in Figure 7c,d, respectively. To match the S n RF definition, an upwelling flux was defined as negative in the present study.  The observed annual mean TOA upward shortwave flux changes (Figure 7b) were highly consistent, both spatially and temporally, with SnRF anomalies (Figure 6b) across the NH during 2000-2019. Similar to SnRF, the TOA upward shortwave flux increased at high latitudes in NA and EU, indicating less TOA shortwave flux back to space and more shortwave radiation absorption by the Earth climate system from 2000 to 2019. In comparison, the TOA upward shortwave flux was enhanced in the mid-latitudes of Central Asia, northeastern China, and Central United States during 2000-2019, which 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 large-scale cooling The observed annual mean TOA upward shortwave flux changes (Figure 7b) were highly consistent, both spatially and temporally, with S n RF anomalies (Figure 6b) across the NH during 2000-2019. Similar to S n RF, the TOA upward shortwave flux increased at high latitudes in NA and EU, indicating less TOA shortwave flux back to space and more shortwave radiation absorption by the Earth climate system from 2000 to 2019. In comparison, the TOA upward shortwave flux was enhanced in the mid-latitudes of Central Asia, northeastern China, and Central United States during 2000-2019, which 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 large-scale cooling trends for land surface temperature in winter over at mid-latitudes [46,47]. Meanwhile, the interannual variability of the annual mean TOA upward shortwave flux increased at a rate of 0.181 Wm −2 decade −1 (p < 0.05), which is consistent with the estimated S n RF anomalies during the same period (Figure 6c). The correlation coefficient between the S n RF and TOA upward shortwave flux in the NH was 0.79 (p < 0.05), highlighting the contribution of S n RF to TOA upward shortwave anomalies in the NH from 2000 to 2019. The strong correlation between S n RF and TOA upward shortwave flux indicates that changes in S n RF may further influence and control the local and regional energy budgets owing to climate change.

Attribution of S n RF Anomalies over the NH during 2000-2019
Changes in the annual mean S n RF could be positively represented by S n RF every month by utilizing a multiple linear regression equation with different coefficients, as displayed in Equation (5). To attribute the annual mean S n RF changes and explore the dominant month that controlled the annual mean S n RF, we calculated the relative contribution of S n RF in each month to the annual mean S n RF anomalies over the NH from 2000 to 2019. Applying a ridge regression analysis with Equation (5), the relative contribution of S n RF in each month to the annual mean S n RF anomalies over the NH during 2000-2019 was calculated, as displayed in Figure 8. Based on the results of SCE and SnRF anomalies over the NH, we inferred that the annual mean SnRF over the NH exhibited a continuous weakening trend from 2000 to 2019, despite the negligible changes in the annual mean SCE. With a rapid SCE decline in snow melting seasons, the SnRF weakened during the melting seasons, especially in May. Conversely, although SCE increased in accumulation seasons, the low insolation and solar incidence angle led to a low SnRF, contributing little to the annual mean SnRF. The Z-score of the annual mean S n RF over the NH during 2000-2019 is displayed in Figure 8a. The resulting contributions from the S n RF in each month to the annual mean S n RF anomalies were significantly different, wherein the S n RF in melting seasons controlled the annual mean S n RF variability from 2000 to 2019. Detailed 20-year averaged relative contributions from each month to the annual mean S n RF anomalies (Figure 8b) further revealed that May accounted for 23.17% of the annual mean S n RF variability over the NH, followed by March (18.11%) and April (16.82%) during 2000-2019. The S n RF in accumulation seasons contributed insignificantly towards the annual mean S n RF from 2000 to 2019. Although the SCE over the NH significantly increased in October (Figure 3a), the S n RF in October remained low because of low insolation during that month.

Discussion
Based on the results of SCE and S n RF anomalies over the NH, we inferred that the annual mean S n RF over the NH exhibited a continuous weakening trend from 2000 to 2019, despite the negligible changes in the annual mean SCE. With a rapid SCE decline in snow melting seasons, the S n RF weakened during the melting seasons, especially in May. Conversely, although SCE increased in accumulation seasons, the low insolation and solar incidence angle led to a low S n RF, contributing little to the annual mean S n RF.

Discussion
Snow-induced radiative forcing S n RF quantitatively estimates the influence of snow on the energy balance of the Earth-atmosphere system and reflects the importance of snow in the mechanism of potential climate change. Therefore, estimating the S n RF is valuable for climate change and climate projections studies. Consequently, quantifying the S n RF and providing support for accurate simulations of the climate change response were the major objectives of the present study.
Previous studies have proved that SCE in the NH has experienced a rapid decrease since satellite observations began in the late 1960s [3,[36][37][38]. Climate projections suggest that the SCE will continue shrinking in the future [37,39], coincident with hemispheric warming, indicating the positive feedback of surface reflectivity on climate [2,36,37]. Differently to the SCE projections, the rapid SCE decline over the NH has experienced a hiatus in recent years, especially during 2000-2019. Therefore, novel investigations into SCE variability and its contribution to the Earth energy budget are necessary. Due to the limited spatial coverage of the surface albedo dataset, high latitudes above 75 • N were excluded in this study.
The method utilized in this study is appropriate for optimizing the input parameters (e.g., snow cover, surface albedo, and albedo radiative kernels) in S n RF calculations to enhance climate change studies and climate projections in the NH. Although we employed the best achievable satellite observations, several issues and limitations exist. Unlike previous research, this study utilized C3S ALBH-DH in S n RF estimation, which provides a better correlation coefficient with SCE anomalies, both in March and June, during 2000-2019. Using C3S ALBH-DH in S n RF estimation greatly reduced the uncertainty in the S n RF calculation induced by surface albedo biases that exist at coarser spatial resolution, such as in GLASS [5] and MODIS [2,7]. To verify the performance of the C3S ALBH-DH climate studies, we utilized GLASS and CLARA-SAL-A2 for a comparison. As CLARA-SAL-A2 solely provides black-sky surface albedo (integration of the bidirectional reflectance distribution function over the viewing hemisphere), the cross-comparison analysis focused on this. Consequently, S n RF in this study only represents the shortwave (solar spectrum) component of the snow-induced radiative forcing. However, compared with S n RF using white-sky surface albedo (integration of the directional albedo over the illumination hemisphere), the comparable value using black-sky surface albedo has a slightly larger value of 4.32% [2,9]. As the difference between black-sky and white-sky surface albedo is systematic, this difference is not likely to influence the spatial and temporal changes in S n RF.
The albedo radiative kernel approach facilitates the estimation of S n RF in this study. However, unsolved spectral variations were excluded from this study. Snow induces a greater albedo contrast in the visible region than that in the near-infrared region, whereas the albedo radiative kernels are generated with spectrally uniform surface albedo perturbations, and with spectrally varying surface fluxes. A sensitivity analysis by Flanner et al. [2] observed that unsolved spectral variations may result in a 5% difference between the results from the entire broadband and separate components. Thus, the excluded unresolved spectral variations in the radiative kernel technique did not influence the major conclusions of our study.
Except for S n RF in clear-sky conditions, the comparable results in all-sky conditions also require further attention. However, published studies using both all-sky and clear-sky albedo radiative kernels have reported that the all-sky S n RF exhibits similarly inter-annual, intra-annual, and spatial pattern variabilities to clear-sky S n RF, except for a 30% difference in magnitude [2,11,48]. Therefore, this study solely utilized clear-sky albedo radiative kernels in S n RF calculation, which is not likely to influence the major conclusions of this study.

Conclusions
Based on satellite observations, this study explored the latest SCE changes over the NH, quantified S n RF, and determined its influence on TOA shortwave anomalies from 2000 to 2019. These results provide insights into the role of snow cover in the Earth's climate system in the context of climate change.
A hiatus in the rapid decline of SCE has been observed in recent years, especially during 2000-2019. In this study, detailed SCE analysis in each month revealed contrasting SCE anomalies between the accumulation and melting seasons over the NH during 2000-2019 compared with those during 1972-2000. A cross-comparison between the snow cover percentage and three surface albedo values in March and June, namely those from C3S ALBH-DH, CLARA-SAL-A2, and GLASS, from 2000 to 2019, demonstrated the reliability of C3S ALBH-DH for estimating S n RF. Compared with GLASS (r = 0.72, p < 0.05) and CLARA-SAL-A2 (r = 0.64, p < 0.05), the C3S ALBH-DH (r = 0.76, p < 0.05) displayed a better linear correlation coefficient with snow cover percentage in March. Meanwhile, the performance of C3S ALBH-DH (0.49, p < 0.05) in June was better than that of GLASS (r = 0.48, p < 0.05) and CLARA-SAL-A2 (r = 0.39, p > 0.05) during 2000-2019.
Utilizing the MOD10CM snow cover range, C3S ALBH-DH surface albedo, and three albedo radiative kernels, this study estimated the S n RF from 2000 to 2019. The 20-year annual mean S n RF in the NH was estimated to be −1.13 (±0.05) W m −2 from 2000 to 2019. Meanwhile, changes in S n RF during 2000−2019 were estimated to be 0.094 (±0.0081) W m −2 , indicating a generally weakening trend of radiative cooling over the NH. The coefficient of determination between S n RF and the TOA upward shortwave flux over the NH from 2000 to 2019 was 0.79 (p < 0.05), highlighting the contribution of S n RF to TOA upward shortwave anomalies over the NH. In addition, the ridge regression analysis revealed that S n RF in the snow accumulation seasons dominated the annual mean S n RF variability over the NH during 2000-2019. May accounted for the greatest (23.17%) annual mean S n RF variability over the NH, followed by March (18.11%) and April (16.82%).
A generalized analysis observed that, although the annual mean SCE exhibited negligible changes from 2000 to 2019, the notable seasonal SCE changes led to a continuous weakening of S n RF over the NH, which should be considered when performing climate projection studies. Therefore, the results presented in this study are beneficial for climate change studies. Owing to the continuous seasonal variability and spatial heterogeneity of SCE in the NH, the monitoring, modeling, and prediction of S n RF are urgently required to reduce their bias in climate projections and CMIP models.