Satellite-based radiative forcing by light-absorbing particles in snow across the Northern Hemisphere

Snow is the most reflective natural surface on Earth and consequently plays an important role in Earth's climate. Light-absorbing particles (LAPs) deposited on the snow surface can effectively decrease snow albedo, resulting in positive radiative forcing. In this study, we used remote-sensing data from NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) and the Snow, Ice, and Aerosol Radiative (SNICAR) model to quantify the reduction in snow albedo due to LAPs before validating and correcting the data against in situ observations. We then incorporated these corrected albedo-reduction data in the Santa Barbara DISORT (Discrete Ordinate Radiative Transfer) Atmospheric Radiative Transfer (SBDART) model to estimate Northern Hemisphere radiative forcing except for midlatitude mountains in December–May for the period 2003–2018. Our analysis reveals an average corrected reduction in snow albedo ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$) of ∼ 0.021 under all-sky conditions, with daily radiative forcing (${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$) values of ∼ 2.9 W m⁻², over land areas with complete or near-complete snow cover and with little or no vegetation above the snow in the Northern Hemisphere. We also observed significant spatial variations in $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$, with the lowest respective values (∼ 0.016 and ∼ 2.6 W m⁻²) occurring in the Arctic and the highest (∼ 0.11 and ∼ 12 W m⁻²) in northeastern China. From MODIS retrievals, we determined that the LAP content of snow accounts for 84 % and 70 % of the spatial variability in albedo reduction and radiative forcing, respectively. We also compared retrieved radiative forcing values with those of earlier studies, including local-scale observations, remote-sensing retrievals, and model-based estimates. Ultimately, estimates of radiative forcing based on satellite-retrieved data are shown to represent true conditions on both regional and global scales.

1 Introduction

Seasonal snow cover affects 30 % of Earth's land surface and exerts a cooling influence on global climate through its direct interaction with the surface radiance budget (Painter et al., 1998; Flanner et al., 2011). However, snow surface darkening due to light-absorbing particles (LAPs) such as black carbon (BC), organic carbon (OC), dust, and algae can significantly alter the reflective properties of snow (Warren, 1982, 1984; Hadley and Kirchstetter, 2012). When deposited on the snow surface, LAPs increase the absorption of solar radiation (Painter et al., 2012a; Liou et al., 2014; Wang et al., 2014, 2019; Dang et al., 2017), thereby reducing the snow albedo (Warren and Brandt, 2008; Kaspari et al., 2014). As a result, radiative forcing of LAPs in snow (RFLS) plays a critical role in snow-cover decline on both regional and global scales (Wiscombe and Warren, 1980), perturbing the climate system and impacting hydrological cycles (Barnett et al., 2005; Qian et al., 2011).

One of the primary LAPs, BC, is derived from the incomplete combustion of fossil fuels and biomass (Bond et al., 2013; Dang et al., 2015) and is second only to CO₂ in its contribution to climate forcing (Hansen and Nazarenko, 2004; Ramanathan and Carmichael, 2008; Bond et al., 2013). Yet, despite considerable efforts to measure the BC content of Northern Hemisphere snow and ice (Doherty et al., 2010, 2014; Huang et al., 2011; Ye et al., 2012; Wang et al., 2013, 2017; Shi et al., 2020), the inherent challenges presented by a temporospatially variable snow cover mean our understanding of LAPs in snow is far from complete. As a result, persistent uncertainties remain in regional- and global-scale RFLS estimates based on field measurements (Zhao et al., 2014).

Several previous investigations have utilized numerical models to estimate RFLS, including that of Hansen and Nazarenko (2004), who concluded that BC in snow and ice exerts a positive climate forcing throughout the Northern Hemisphere of +0.3 W m⁻², explaining approximately one-fourth of observed global warming. More recently, Flanner et al. (2007) employed an aerosol–chemical-transport general-circulation model, coupled with the Snow, Ice, and Aerosol Radiative (SNICAR) model (Flanner et al., 2007, 2009), to estimate globally averaged radiative forcing values of +0.054 (range of 0.007–0.13) and +0.049 (0.007–0.12) W m⁻² for a strong (1998) and weak (2001) boreal fire years, respectively. Using the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008) coupled with a chemistry component (Chem) (Grell et al., 2005) and SNICAR modeling, Zhao et al. (2014) demonstrated that RFLS over northern China in January–February 2010 was ∼ 10 W m⁻². However, despite their potentially valuable contribution, climate models contain significant uncertainties in representations of LAP emissions, transport, deposition, and post-depositional processes that can propagate into simulations of LAP concentrations and their climate forcing (Qian et al., 2015; Lee et al., 2016). Zhao et al. (2014) also confirmed that, relative to observational data, modeled LAPs and radiative forcing estimates exhibit biases that are difficult to explain and quantify. These shortcomings underscore the need for a refined approach to estimating real-time RFLS that minimizes the mismatch between field observations and model simulations.

In addition to modeling, remote sensing has been used to assess the physical characteristics of snow cover (Nolin and Dozier, 1993, 2000; Painter et al., 2009, 2012a, 2013; Miller et al., 2016). Nolin and Dozier (2000), for example, retrieved grain size data from satellite-derived reflectance at near-infrared (NIR) wavelengths following the rationale that snow-grain size, in conjunction with solar zenith angle, dictates the path length of penetrating photons (Wiscombe and Warren, 1980) and thus influences albedo in the NIR. Similarly, recent studies have attempted to employ satellite-derived snow albedo at visible (VIS) wavelengths to retrieve RFLS data (Seidel et al., 2016; Pu et al., 2019). Briefly, this retrieval method exploits the imaginary component of the complex refractive index for ice (K_ice), which is very low at VIS wavelengths and results in the extremely high VIS albedo for pure snow. In contrast, the imaginary component of the complex refractive index for LAPs (K_LAPs) at VIS wavelengths is orders of magnitude greater, resulting in the reduction in VIS snow albedo (Wiscombe and Warren, 1980). Moreover, albedo variability at VIS wavelengths is dominated by even minor concentrations of LAPs (Brandt et al., 2011; Painter et al., 2012b).

Painter et al. (2012a) employed surface-reflectance data provided by NASA's Moderate Resolution Imaging Spectroradiometer (MODIS) for the upper Colorado River basin and Hindu Kush Himalaya (HKH) to make the first quantitative, remote-sensing-based retrievals of instantaneous surface radiative forcing (RF) due to LAPs. Relative to the Western Energy Balance of Snow (WEBS) network (Painter et al., 2007), that study established that MODIS-derived radiative forcing exhibits a positive bias at lower RF values and a slightly negative bias at higher values. A more recent study by Seidel et al. (2016) used remote sensing to constrain instantaneous melt-season RFLS values of 20–200 W m⁻² for the Sierra Nevada and Rocky Mountains, while Pu et al. (2019) reported MODIS-derived values of 22–65 W m⁻² for northern China in January–February (regional average of ∼ 45 W m⁻²). Acknowledging this demonstrated efficacy of remote-sensing retrievals for establishing RFLS on regional scales, we note that this approach has so far not captured spatial variability in RFLS on a global scale.

In this study, we employed MODIS data to determine the reduction in the Northern Hemisphere snow albedo due to LAPs. Retrievals were validated and corrected according to ground-based snow observations, after which spatial variability in albedo reduction and radiative forcing over mapped snow-covered area in the Northern Hemisphere was assessed quantitatively. Finally, we compared our satellite-derived radiative forcing values with the modeling results of CESM2 (Eyring et al., 2016; Danabasoglu et al., 2020). Despite the persistence of non-negligible uncertainties and biases, our satellite-based retrievals constitute the first hemisphere-scale assessment of RFLS and provide valuable information for improving climate model simulations.

2 Data

2.1 Remote-sensing data

To investigate the impact of LAPs on snow albedo, we utilized the following MODIS data sets: surface albedo (MCD43C3; 0.05^∘ × 0.05^∘ resolution), snow cover (MYD10C1; 0.05^∘ × 0.05^∘ resolution), land cover type (MCD12C1; 0.05^∘ × 0.05^∘ resolution), and atmospheric parameters (MYD08_D3; 1^∘ × 1^∘ resolution). Each data set corresponds to December–May for the period 2003–2018 (https://earthdata.nasa.gov, last access: 20 January 2019). MCD43C3 is the daily combined MODIS output derived from both the Terra and Aqua satellites and provides black-sky albedo (directional hemispherical reflectance, DHF) and white-sky albedo (bi-hemispherical reflectance, BHF) at local solar noon for bands 1–7 (band 1, 620–670 nm; band 2, 841–876 nm; band 3, 459–479 nm; band 4, 545–565 nm; band 5, 1230–1250 nm; band 6, 1628–1652 nm; band 7, 2105–2155 nm) as well as values for quality control, local noon solar zenith angle, and associated parameters. MCD43C3 observations are weighted to estimate albedo on the ninth day of each 16 d period and have been corrected for the influence of local slope and aspect, atmospheric gases, and aerosols.

Snow-cover data are provided daily by MYD10C1 as a report of the snow-cover fraction (SCF), derived from the Normalized Difference Snow Index (NDSI). MCD12C1 provides a spatially aggregated and reprojected land cover type, which is derived from the supervised classification of MODIS reflectance data, while MODIS MYD08_D3 reports values of solar azimuth angle.

Daily average solar radiances and cloud fraction were obtained from NASA's Clouds and the Earth's Radiant Energy System (CERES: https://ceres.larc.nasa.gov, last access: 12 April 2019), part of the Earth Observing System comprising the Aqua, Terra, and S-NPP satellites. CERES provides instantaneous measurements of solar radiances, which are then converted to daily average flux by angular dependence and empirical diurnal albedo modeling as the satellite passes through the point of descent (Doelling et al., 2013; Su et al., 2015; Loeb et al., 2018). We used the total downward shortwave flux and cloud fraction at the surface, provided by the “CERES Single Scanner Footprint 1.08 (SSF1deg)” product, to estimate daily average RFLS under all-sky conditions.

Shuttle Radar Topography Mission (SRTM) digital elevation data are provided by the US Geological Survey (https://www.usgs.gov/, last access: 9 December 2018) to adjust slope- and aspect-induced changes in surface solar irradiance in complex terrain. The spatial resolution of SRTM data for the Northern Hemisphere is 30 m.

2.2 Snow depth data

Estimates of snow depth were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Reanalysis (ERA-Interim) (https://www.ecmwf.int, last access: 15 January 2019). ERA-Interim is a new generation of reanalysis based on a 12-hourly and four-dimensional variational data assimilation (4D-Var) covering the period 1979–present. ERA-Interim performs better in model physics frameworks, data quality control, and background error criteria than previous versions (Berrisford et al., 2011; Brun et al., 2013). In this study, we used snow-water equivalent (SWE) data for December–May covering the period 2003–2018. These data were generated by forecast models and updated according to a Cressman analysis of snow observations (Drusch et al., 2004; Dee et al., 2011). We note that the previous occurrence of false snow-free patches, arising from the application of Cressman analysis in regions of sparse ground control, has been mitigated by ECMWF upgrades (Dee et al., 2011). Finally, SWE is converted to snow depth by assuming that average December–May snow density is ∼ 300 kg m⁻³, consistent with snow depth estimates by the Canadian Meteorological Centre (CMC) (Sturm et al., 1995; Brown and Mote, 2009).

2.3 In situ measurements of LAPs in snow

To correct the satellite retrievals, we collected a comprehensive set of in situ measurements of BC concentrations from the field campaigns in the Arctic in spring of 2005–2009 (Doherty et al., 2010); North America in January–March of 2013 (Doherty et al., 2014); northern China in January–February of 2010, 2012, and 2014 (Ye et al., 2012; Wang et al., 2013, 2017). The BC concentrations are measured by the two-sphere integrating-sandwich (TSI) spectrophotometer in the Arctic, North America, and northern China (Grenfell et al., 2011; Wang et al., 2020). Briefly, TSI produces a diffuse radiation field when the white light illumination is transmitted into an integrating sphere; then the diffuse radiation passes through the filter and is detected by a spectrometer. The TSI technique acquires the light attenuation spectrum due to the LAPs loaded on the sample filter (Grenfell et al., 2011). Then, the light attenuation spectrum of the sample filter is transformed into an equivalent BC mass (unit: g cm⁻²) loading by comparing against the standard filters. The equivalent BC has been defined by Doherty et al. (2010), which briefly as the amount of BC in the snow accounted for the wavelength-integrated total light absorption in the wavelengths of 300–750 nm by all particulate constituents. In this study, we used BC_equiv for all LAPs to calculate the in situ snow albedo reduction and radiative forcing (Fig. S3).

2.4 Climate model simulations

We compared our remotely sensed retrievals of daily average RFLS for the 2003–2014 study period with simulated results derived from CESM2 (https://esgf-node.llnl.gov/, last access: 15 July 2019). In this study, we employed simulations of snow BC concentrations derived from the CESM2 historical experiments in conjunction with ERA-Interim SWE, MODIS-retrieved snow grain size, and CERES total downward shortwave flux data under all-sky conditions, to model daily average RFLS for the study period. Simulations were performed using the Snow, Ice, and Aerosol Radiative (SNICAR) and Santa Barbara DISORT (Discrete Ordinate Radiative Transfer) Atmospheric Radiative Transfer (SBDART) models, and the model output was compared with satellite-based retrievals.

3 Methods

3.1 Radiative-transfer model

In this study, we used the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) model to calculate spectral surface solar irradiance. Constituting one of the most widely applied models for calculating the atmospheric radiative transfer at Earth's surface, under both clear- and cloudy-sky conditions (Ricchiazzi et al., 1998), SBDART combines a low-resolution atmospheric transmission model, Discrete Ordinate Radiative Transfer (DISORT) module, and Mie scattering output for the scattering of light by ice crystals and water droplets (Stamnes et al., 1988; Fu et al., 2017). Radiative-transfer equations for a vertically inhomogeneous, non-isothermal, plane-parallel atmosphere are integrated numerically using the DISORT module. SBDART comprises multiple standard atmospheric profiles, cloud models, and basic surface types as well as vertical-distribution models for aerosols and gas absorption and enables users to specify these input parameters in real values. In our study, the subarctic and midlatitude winter standard atmospheric conditions are assumed, and the tropospheric and stratospheric background aerosols are archived in SBDART (Tanre et al., 1990). According to Dang et al. (2017), the cloud optical depth at high latitudes and midlatitudes was assumed to be 11 and 20 under cloudy-sky conditions, respectively. The spectral irradiance from SBDART is only used for integrating the spectral MODIS albedo to achieve broadband albedo; thus the uncertainty in solar irradiance from the assumed atmospheric properties has limited influence on the retrieval of radiative forcing (see Sect. 3.2). Average incident direct and diffuse solar spectra for December–May under clear and cloudy skies are shown in Fig. S1.

The Snow, Ice, and Aerosol Radiative (SNICAR) model is a two-stream multiple-scattering radiative-transfer model (Flanner et al., 2007, 2009) that has been used widely both to simulate the albedo, transmission, and vertical absorptivity of LAP-contaminated snowpack and to estimate RFLS (Painter et al., 2012a; Bryant et al., 2013; Miller et al., 2016). SNICAR employs the theory proposed by Wiscombe and Warren (1980) and Toon et al. (1989). Specifically, snow is considered to be composed of aggregated ice grains with optical effective radii (R_eff) of 50–1500 µm, lognormal distribution, and spherical grain shape. SNICAR also accounts for the incident radiation at the surface and its spectral distribution, solar zenith angle, snow depth and density, snow layer number, and the type and concentration of LAPs in the snowpack. The model's ability to provide realistic simulations of snow albedo has been verified by several previous studies (Hadley and Kirchstetter, 2012; Meinander et al., 2013; Zhong et al., 2017; Wang et al., 2017).

3.2 Retrieval of quantitative snow properties from remote sensing

The variability in spectral snow albedo depends on the LAP content, grain size, grain shape, and depth of the snowpack, in addition to solar zenith angle. As shown in Fig. 1a, the deposition of BC (as representative of LAPs generally) serves to decrease the albedo of snow significantly, particularly in the ultraviolet (UV) and VIS wavelengths, which account for approximately half of all direct solar irradiance and the majority of diffuse solar irradiance (Fig. S1). In contrast, the impact of BC on albedo is considerably smaller in NIR wavelengths and can be negligible at > ∼ 1150 nm. Snow depth plays a similar role to LAP content and primarily affects albedo in UV and VIS wavelengths (Fig. 1b).

Figure 1 Variations in spectral snow albedo due to (a) LAP content (ng g⁻¹), (b) snow depth (m), (c) snow grain size (µm), and (d) solar zenith angle (^∘) while the other three parameters are kept constant.

Although snow albedo decreases with snow depth, previous studies have tended to assume a semi-infinite snowpack for which albedo is independent of depth. As a consequence, the role of LAPs in albedo reduction has been overestimated for those areas where the snowpack is thin (Warren, 2013). In this study, we incorporated ERA-Interim SWE data in our SNICAR model simulations to correct for the snow depth overestimation effect. In contrast, snow grain size and solar zenith angle influence the snow albedo chiefly in NIR wavelengths (Fig. 1c, d). Specifically, albedo tends to decrease with increasing snow grain size and declining solar zenith angle. In this study, we derived quantitative snow parameters (grain size, albedo reduction, and RFLS) from MODIS data in conjunction with the SNICAR and SBDART models. The specific workflow for retrieving RFLS from satellite data is shown in Fig. 2.

Figure 2 Workflow depicting the calculation and validation of radiative forcing of LAPs in snow: the blue boxes denote the external input data, while the orange boxes are used for calculations in this study.

3.2.1 Retrieval of blue-sky albedo

MCD43 provides black-sky and white-sky albedo, which are defined as albedo in the absence of the diffuse and direct competent of solar irradiance. Accordingly, the actual spectral albedo for a land surface at wavelength λ (also called blue-sky albedo: ${\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{blue-clear}}$) under clear-sky conditions can be calculated as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{blue-clear}}=f_{\mathrm{dif},\mathit{\lambda }}^{\text{clear}}\times {\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{white-sky}}+(\mathrm{1}-f_{\mathrm{dif},\mathit{\lambda }}^{\text{clear}})\times {\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{black-sky}},\end{array}$$

where ${\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{white-sky}}$ and ${\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{black-sky}}$ are MODIS-derived values for white-sky and black-sky albedo, respectively, and $f_{\mathrm{dif},\mathit{\lambda }}^{\text{clear}}$ is the ratio of diffuse irradiance to the total solar irradiance under clear skies (Lewis and Barnsley, 1994). The latter is calculated as follows:

$$\begin{array}{}\text{(2)}& f_{\mathrm{dif},\mathit{\lambda }}^{\text{clear}}={\displaystyle \frac{E_{\mathrm{dif}}^{\text{clear}}(\mathit{\lambda };\mathit{\phi })}{E_{\mathrm{dif}}^{\text{clear}}\left(\mathit{\lambda };\mathit{\phi }\right)+E_{\mathrm{dir}}^{\text{clear}}(\mathit{\lambda };\mathit{\phi })\times \cos \mathit{\beta }}},\end{array}$$

where φ is latitude; $E_{\mathrm{dif}}^{\text{clear}}\left(\mathit{\lambda };\mathit{\phi }\right)$ denotes the diffuse spectral irradiance on a horizontal surface; and $E_{\mathrm{dir}}^{\text{clear}}(\mathit{\lambda };\mathit{\phi })$ denotes the direct spectral irradiance on a surface perpendicular to the sun, derived from the SBDART model under clear-sky conditions. β represents local solar zenith angle, which is obtained using the topographic correction method (Teillet et al., 1982; Negi and Kokhanovsky, 2011):

$$\begin{array}{}\text{(3)}& \cos \mathit{\beta }=\cos {\mathit{\theta }}_{\mathrm{0}}\cos {\mathit{\theta }}_T+\sin {\mathit{\theta }}_{\mathrm{0}}\sin {\mathit{\theta }}_T\cos ({\mathit{\varphi }}_{\mathrm{0}}-{\mathit{\varphi }}_T),\end{array}$$

for which θ₀ represents the solar zenith angle for a horizontal surface; ϕ₀ is the solar azimuth angle; and θ_T and ϕ_T denote slope inclination and aspect, respectively. Similarly, we can derive the blue-sky albedo for cloudy-sky conditions (${\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{blue-cloudy}}$). Then, we used cloud fraction (f_cloud) from CERES to weight clear-sky albedo and cloudy-sky albedo to obtain actual all-sky albedo (${\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\mathrm{all}}$):

$$\begin{array}{}\text{(4)}& {\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\mathrm{all}}=f_{\text{cloud}}\times {\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{blue-cloudy}}+(\mathrm{1}-f_{\text{cloud}})\times {\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\text{blue-clear}}.\end{array}$$

3.2.2 Retrieval of snow cover and albedo values

As shown in Fig. 2, the snow-covered area is mapped according to the actual all-sky albedo (${\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\mathrm{all}}$) in band 4 (band center ∼ 555 nm) and the Normalized Difference Snow Index (NDSI), both of which are required to exceed 0.6 (Negi and Kokhanovsky, 2011). According to the MODIS Snow Products Collection 6 user guide (http://nsidc.org/data, last access: 1 October 2020), the fractional snow cover (FSC) can be calculated as follows:

$$\begin{array}{}\text{(5)}& \mathrm{FSC}=-\mathrm{0.01}+\mathrm{1.45}\times \mathrm{NDSI}.\end{array}$$

Accordingly, the identified snow-covered area (ISCA) has an FSC value of > 86 % but not always 100 %. Therefore, the MODIS-derived albedo for a particular ISCA is a combination of values representing both snow and the snow-free underlying surface. Following Pu et al. (2019), the snow albedo (${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{all}}$) can be distinguished from the mixed albedo by the equation

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}& {\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\mathrm{all}}=\\ & {\displaystyle \frac{E_{\text{all-sky},\mathit{\lambda }}\times \mathrm{FSC}\times {\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{all}}+E_{\text{all-sky},\mathit{\lambda }}\times \left(\mathrm{1}-\mathrm{FSC}\right)\times {\mathit{\alpha }}_{\mathrm{underlying},\mathit{\lambda }}}{E_{\text{all-sky},\mathit{\lambda }}}}\\ & =\mathrm{FSC}\times {\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{all}}+\left(\mathrm{1}-\mathrm{FSC}\right)\times {\mathit{\alpha }}_{\mathrm{underlying},\mathit{\lambda }}\end{array}\text{(7)}& {\displaystyle }{\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{all}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{MODIS},\mathit{\lambda }}^{\mathrm{all}}-\left(\mathrm{1}-\mathrm{FSC}\right)\times {\mathit{\alpha }}_{\mathrm{underlying},\mathit{\lambda }}}{\mathrm{FSC}}},\end{array}$$

where E_all-sky,λ is total solar irradiance under all-sky conditions, a linear combination of direct and diffuse components of solar irradiance under clear skies and cloudy skies using a similar strategy via Eqs. (1)–(4); α_{underlying,λ} represents the albedo of the underlying surface and was obtained from Siegmund and Menz (2005). As depicted in Fig. 3b, vegetation and bare soil are the main types of underlying surface in the ISCA.

Figure 3 Spatial distributions of (a) the identified snow-covered area (ISCA) and (b) the different land cover types based on MODIS data for the Northern Hemisphere. ISCAs (white) can be separated into northeastern China (NEC), Eurasia (EUA), North America (NA), and the Arctic.

3.2.3 Retrieval of snow grain size

The snow optical-equivalent grain size (R_eff) is retrieved by fitting SNICAR-simulated snow albedo to MODIS-derived snow albedo at 1240 nm (the central wavelength of MODIS band 5) following the protocol of Nolin and Dozier (2000). This retrieval method is not influenced by liquid water and water vapor and has been employed widely in previous studies (e.g., Painter et al., 2013; Seidel et al., 2016). Both Nolin and Dozier (2000) and Pu et al. (2019) reported that the retrieved R_eff compares favorably with ground-based measurements of snow grain size. In this study, we chose to exclude the ISCA, where MODIS-derived snow albedo at 1240 nm is < 0.3, to avoid misrepresenting R_eff (Tedesco and Kokhanovsky, 2007).

3.2.4 Retrieval of snow albedo reduction and RFLS

The spectrally integrated reduction in snow albedo due to LAPs ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,noon}}^{\mathrm{LAPs}}$) is estimated for local-noon and all-sky conditions using solar irradiance and the difference between MODIS-derived spectral snow albedo (${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{all}}$) and simulated pure snow albedo (${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{mdl}}$). Because MODIS provides only four VIS bands, we fitted snow albedo data obtained via MODIS to a continuous 300–2500 nm spectrum (${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{MODIS}}$, with a 10 nm interval) following the method provided by Pu et al. (2019). Thereafter, the broadband albedo reduction due to LAPs retrieved from MODIS ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,noon}}^{\mathrm{LAPs}}$) can be calculated as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,noon}}^{\mathrm{LAPs}}=\\ & {\displaystyle \frac{\sum _{\mathit{\lambda }=\mathrm{300}\mathrm{nm}}^{\mathit{\lambda }=\mathrm{2500}\mathrm{nm}}\left({\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{mdl}}-{\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{MODIS}}\right)\cdot E_{\text{all-sky},\mathit{\lambda }}\cdot \mathrm{\Delta }\mathit{\lambda }}{\sum _{\mathit{\lambda }=\mathrm{300}\mathrm{nm}}^{\mathit{\lambda }=\mathrm{2500}\mathrm{nm}}{E}_{\text{all-sky},\mathit{\lambda }}\cdot \mathrm{\Delta }\mathit{\lambda }}}\end{array},\end{array}$$

where ${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{mdl}}$ is the pure snow albedo simulated by SNICAR using MODIS-derived R_eff and ERA-Interim snow depth data, ${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{MODIS}}$ is the continuous snow albedo derived from MODIS retrievals, and Δλ is 10 nm.

Following Miller et al. (2016), we assumed that the properties for snow and LAPs remain invariable throughout the day. Based on calculated ${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{mdl}}$ and ${\mathit{\alpha }}_{\mathrm{snow},\mathit{\lambda }}^{\mathrm{MODIS}}$ at local noon, the diurnal variation in pure and polluted snow albedo can be simulated by SNICAR from sunrise to sunset. Then, daily average snow albedo reduction ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$) can be derived by integrating the diurnal snow albedo reduction, which is weighted by simultaneous solar irradiance from SBDART. Similarly, we used measurements of LAPs in contaminated snow to calculate the in situ reduction in snow albedo ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{in situ,daily}}^{\mathrm{LAPs}}$). To derive a correction factor for MODIS retrievals, we applied a similar validation strategy to that of Zhu et al. (2017):

$$\begin{array}{}\text{(9)}& c={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left({\displaystyle \frac{\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}}{\mathrm{\Delta }{\mathit{\alpha }}_{\text{in situ,daily}}^{\mathrm{LAPs}}}}\right),\end{array}$$

where c is the correction factor for $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$, and n is the number of the respective in situ measurements. Accordingly, the corrected albedo reduction ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$) is calculated as follows:

$$\begin{array}{}\text{(10)}& \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}={\displaystyle \frac{\mathrm{1}}{c}}\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}.\end{array}$$

The daily average, spectrally integrated RFLS (${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$) is calculated for all-sky conditions as follows:

$$\begin{array}{}\text{(11)}& {\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}=\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\cdot {\mathrm{SW}}_{\text{all-sky}},\end{array}$$

where SW_all-sky represents the daily average total downward shortwave fluxes, obtained from CERES under all-sky conditions.

3.2.5 Attribution of spatial variability in snow albedo reductions and radiative forcing

As demonstrated above, reductions in snow albedo and RFLS are dependent primarily on LAP content, R_eff, snow depth (SD), solar zenith angle, surface topography, and solar irradiance, the latter three of which can be categorized as the geographic factor (G). We used an impurity index (I_LAPs) to represent the LAP content of the snowpack (Di Mauro et al., 2015; Pu et al., 2019) following the equation

$$\begin{array}{}\text{(12)}& I_{\mathrm{LAPs}}={\displaystyle \frac{\ln \left({\mathit{\alpha }}_{\text{snow,band4}}^{\mathrm{all}}\right)}{\ln \left({\mathit{\alpha }}_{\text{snow,band5}}^{\mathrm{all}}\right)}},\end{array}$$

where ${\mathit{\alpha }}_{\text{snow,band4}}^{\mathrm{all}}$ and ${\mathit{\alpha }}_{\text{snow,band5}}^{\mathrm{all}}$ are the MODIS-derived snow albedo values for bands 4 and 5, respectively. We then calculated $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ as follows:

$$\begin{array}{}\text{(13)}& \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}=f(I_{\mathrm{LAPs}}{R}_{\mathrm{eff}},\mathrm{SD},G).\end{array}$$

The spatial variability in snow albedo reduction due to I_LAPs can be expressed as

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)=f(I_{\mathrm{LAPs}}\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}},\stackrel{\mathrm{‾}}{\mathrm{SD}},\stackrel{\mathrm{‾}}{G}),\end{array}$$

where $\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}},\stackrel{\mathrm{‾}}{\mathrm{SD}}$, and $\stackrel{\mathrm{‾}}{G}$ indicate spatial-mean values of R_eff, SD, and G, with $\stackrel{\mathrm{‾}}{G}$ requiring spatially constant values for the solar zenith angle, surface topography, and solar irradiance parameters. The following three equations were applied in a similar manner:

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)=f(\stackrel{\mathrm{‾}}{I_{\mathrm{LAPs}}}{R}_{\mathrm{eff}}\stackrel{\mathrm{‾}}{\mathrm{SD}},\stackrel{\mathrm{‾}}{G})\text{(16)}& {\displaystyle }\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(\mathrm{SD}\right)=f(\stackrel{\mathrm{‾}}{I_{\mathrm{LAPs}}}\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}},\mathrm{SD},\stackrel{\mathrm{‾}}{G})\text{(17)}& {\displaystyle }\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(G\right)=f(\stackrel{\mathrm{‾}}{I_{\mathrm{LAPs}}}\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}},\stackrel{\mathrm{‾}}{\mathrm{SD}},G).\end{array}$$

We then fitted $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ through multiple linear regression:

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}=a\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)\\ & +b\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)\\ & +c\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(\mathrm{SD}\right)\\ & +d\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(G\right),\end{array}\end{array}$$

where $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}$ is the fitted snow albedo reduction, and a, b, c, and d denote the regression coefficients. Figure S3a illustrates how $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}$ can explain 99 % of the variance in $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$. Therefore, the attribution of spatial variance in $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ can be replaced with $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}$, enabling Eq. (18) to be written as follows:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}-\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}}=\\ & a\times \left(\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)-\right.\\ & \left.\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)}\right)\\ & +b\times \left(\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)\right.\\ & \left.-\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)}\right)\\ & +c\times \left(\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(\mathrm{SD}\right)\right.\\ & \left.-\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(\mathrm{SD}\right)}\right)\\ & +d\times \left(\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(G\right)\right.\\ & \left.-\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}\left(G\right)}\right),\end{array}\end{array}$$

where $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}-\stackrel{\mathrm{‾}}{\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS}}^{\text{LAPs,fit}}}$ is the snow albedo reduction anomaly ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,anomaly}}^{\text{LAPs,fit}}$). Then, Eq. (19) can be written as

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,anomaly}}^{\text{LAPs,fit}}=a\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)\\ & +b\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)\\ & +c\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}\left(\mathrm{SD}\right)\\ & +d\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}\left(G\right).\end{array}\end{array}$$

According to Huang and Yi (1991) and Pu et al. (2019), the fractional contribution of LAP content to the variability in snow albedo reduction ($R_{\mathrm{\Delta }\mathit{\alpha }}^{\mathrm{LAPs}}$) can be calculated as

$$\begin{array}{}\text{(21)}& {\displaystyle }\begin{array}{rl}& R_{\mathrm{\Delta }\mathit{\alpha }}^{\mathrm{LAPs}}=\\ & {\displaystyle \frac{\mathrm{1}}{m}}\sum _{j=\mathrm{1}}^m{\displaystyle \frac{{\left(a\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}{\left(I_{\mathrm{LAPs}}\right)}_j\right)}^{\mathrm{2}}}{K_j}}\end{array}\text{(22)}& {\displaystyle }\begin{array}{rl}& K_j={\left(a\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}{\left(I_{\mathrm{LAPs}}\right)}_j\right)}^{\mathrm{2}}\\ & +{\left(b\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}{\left(R_{\mathrm{eff}}\right)}_j\right)}^{\mathrm{2}}\\ & +{\left(c\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}{\left(\mathrm{SD}\right)}_j\right)}^{\mathrm{2}}\\ & +{\left(d\times \mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected,anomaly}}^{\mathrm{LAPs}}{\left(G\right)}_j\right)}^{\mathrm{2}},\end{array}\end{array}$$

where m denotes the length of the data set. Values for $R_{\mathrm{\Delta }\mathit{\alpha }}^{R_{\mathrm{eff}}}$, $R_{\mathrm{\Delta }\mathit{\alpha }}^{\mathrm{SD}}$, and $R_{\mathrm{\Delta }\mathit{\alpha }}^G$ can be derived in the same way. Similarly, we can obtain the fractional contribution for daily radiative forcing ($R_{\mathrm{RF}}^{\mathrm{LAPs}}$, $R_{\mathrm{RF}}^{R_{\mathrm{eff}}}$, $R_{\mathrm{RF}}^{\mathrm{SD}}$, and $R_{\mathrm{RF}}^G$).

4 Results

4.1 Study area

Figure 3a depicts the ISCA employed in this study. Most are located in Eurasia, North America, and the Arctic, which are dominated by grassland, shrublands, and bare-soil surfaces (Fig. 3b). Several mid–high-latitude regions that typically support a deep snowpack, including southern Russia, western Europe, and the eastern US, are not identified by MODIS as ISCA due to the broad distributions of forest in those areas (Fig. 3b). This pattern is supported by Bond et al. (2006), who demonstrated that, under such vegetated conditions, LAPs in snow exert a relatively minor influence on radiative forcing. On the other hand, the snowpack over midlatitude mountains at such a coarse resolution (0.05^∘ × 0.05^∘) is too low to identify. In addition, midlatitude mountains are characterized as complex terrain, which will lead to high biases in radiative forcing retrieval at the coarse resolution in spite of topographic correction. Therefore, we did not report the results over midlatitude mountains in this study.

As illustrated in Fig. 3a, ISCAs can be separated into four general regions according to geographical distribution and pollution conditions (Fig. S2a, b): northeastern China (NEC), Eurasia (EUA), North America (NA), and the Arctic. The following analysis of snow albedo reduction and RFLS only concerns ISCA, and the results mainly represent winter for midlatitudes (because spring is mostly snow-free) and spring for the Arctic (because albedos cannot be derived during polar night).

4.2 Global characteristics

Previous studies have highlighted the dominant role of BC in light absorption by snow (Wang et al., 2013; Dang et al., 2017). The spatial distribution of BC emissions density for the Northern Hemisphere in December–May is shown in Fig. S2a. Emissions density exhibits a strong spatial inhomogeneity, ranging from < 10⁻¹ to > 10⁴ g km⁻² per month over ISCA. The highest values occur in NEC, where the emissions are considerably higher than EUA and NA, and the lowest values occur in the Arctic. The wet and dry deposition of BC constitute the primary mechanisms for BC accumulation in snow. As shown in Fig. S2b, the distribution of BC deposition (i.e., the sum of dry and wet deposition) is similar to BC emissions density, with the highest and lowest regional averages corresponding to NEC and the Arctic, respectively. Together, these data indicate that the NEC snowpack is heavily polluted, and thus snow albedo reduction is likely to be highest, while the Arctic snowpack is the least contaminated.

In addition to LAP content, the physical properties of the snowpack, such as depth and grain size, also impact snow albedo (Fig. 1). As depicted in Fig. 4a, the average snowpack in EUA (0.15 m thick) is thinner than in both NA (0.24 m) and NEC (0.19 m), implying a greater impact of snow depth on snow albedo and radiative forcing in EUA. The greatest snow depths occur in the Arctic (> 1 m) and can be considered semi-infinite, meaning that the impact of depth on albedo and radiative forcing is negligible. Figure 4b shows the spatial distribution of MODIS-derived snow grain radius (R_eff). In contrast to BC emissions density, BC deposition, and snow depth, R_eff exhibits minor spatial variability, with regional-average values for NEC, EUA, NA, and the Arctic of 237, 227, 237, and 215 µm, respectively. These values align with the findings of several previous studies (Painter et al., 2013; Seidel et al., 2016; Pu et al., 2019) and imply that the contribution of R_eff to spatial variability in snow albedo reduction and radiative forcing is negligible.

Figure 4 Spatial distributions of 2003–2018 averaged (a) snow depth from ERA-Interim, (d) snow grain size retrieved by MODIS, and (c) total downward shortwave flux at the surface during December–May from CERES.

According to Eq. (11), local solar radiation is an important factor for determining RFLS. Figure 4c depicts the December–May averaged total downward surface shortwave flux under all-sky conditions. Average solar radiative flux values for EUA and NA are comparable to one another but high relative to NEC, which lies at a generally higher latitude (> 40^∘). The lowest values occur in the Arctic due to that region's extreme latitude. The Arctic goes through the polar night during winter so that the radiative effect of LAPs in the Arctic mainly appears in spring. Figure S2d shows the March–May averaged downward surface shortwave flux. As can be seen, the values in the Arctic in March–May are higher than those at midlatitudes in December–February (Fig. S2c). We note that snow albedo reduction and radiative forcing are only calculated over the period when snow-covered area was mapped, which implies that the RFLS will be higher in the Arctic than at midlatitudes for the same snow albedo reduction.

4.3 Corrections based on in situ observations

Albedo reduction calculated using in situ observed LAPs ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{in situ,daily}}^{\mathrm{LAPs}}$) was used to quantitatively correct MODIS retrievals through comparison with MODIS-retrieved snow albedo reduction ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$). Figure S4 displays scatterplots of the ratios of $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ to $\mathrm{\Delta }{\mathit{\alpha }}_{\text{in situ,daily}}^{\mathrm{LAPs}}$ ($r_{\text{in situ}}^{\text{MODIS}}$) for each sampling site (Ye et al., 2012; Wang et al., 2013, 2017; Doherty et al., 2010, 2014). Briefly, for NA, EUA, and the Arctic, where the snowpack is relatively clean, the values for $r_{\text{in situ}}^{\text{MODIS}}$ mostly range between 2 and 10. In contrast, the heavily polluted snowpack in NEC returns $r_{\text{in situ}}^{\text{MODIS}}$ values ranging from 0.5 to 2.5, indicating a negative correlation between the biases of $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ and snow contamination, thus supporting the findings of previous studies (Painter et al., 2012a; Pu et al., 2019). To improve the quality of MODIS retrievals, we developed the correction factors for different regions. According to Eq. (10), the correction factors for NEC, EUA, NA, the Canadian Arctic, the Russian Arctic, and Greenland are 1.6, 4.1, 4.1, 4.4, 5.4, and 6.0, respectively. Hereafter, our analyses are based on the corrected MODIS retrievals. Figure 5 compares the corrected MODIS retrievals to measurement-based results, and the mean absolute error (MAE), root mean square error (RMSE), and correlation coefficient of $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ relative to $\mathrm{\Delta }{\mathit{\alpha }}_{\text{in situ,daily}}^{\mathrm{LAPs}}$ are given in Table S1. Together, these results imply that the corrected MODIS retrievals are plausible. Nevertheless, we note that the correction used in this study is spatially rough due to the low density of in situ measurements and thus that both the uncertainty and bias are non-negligible. To address this issue, we presented further discussion about the accuracy of radiative forcing retrievals (see Sect. 4.5). We also conducted a comprehensive series of comparisons between the MODIS-derived retrievals and values provided via surface measurements, model simulations, and remote sensing (see Sect. 5). We concluded that further field-based measurements of snow albedo are required to improve the quality of satellite retrievals.

Figure 5 Scatterplots of $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ versus $\mathrm{\Delta }{\mathit{\alpha }}_{\text{in situ,daily}}^{\mathrm{LAPs}}$. Panels (a)–(f) represent the snow samples collected in the Canadian Arctic, Russian Arctic, Greenland, North America, northwestern China, and northeastern China, respectively.

4.4 Spatial distributions of snow albedo reduction and radiative forcing

Figure 6a shows the spatial distributions of MODIS-based albedo reduction and daily radiative forcing, and statistics are shown in Fig. 6b and Table 1. On average, $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ provide respective values of 0.021 and 2.9 W m⁻² for Northern Hemisphere ISCA. The highest $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ occurs in NEC, where the regional average of ∼ 0.11 exceeds those of EUA (∼ 0.031) and NA (∼ 0.027) by a factor of ∼ 3–4. This feature reflects the relatively high rate of emissions over NEC, which results in the highest level of BC deposition over ISCA (Fig. S2a, b). In contrast, being located far from major sources of pollution, the relatively clean Arctic snowpack returns the lowest $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ (∼ 0.016) of the entire Northern Hemisphere. Consistent with snow albedo reduction, the highest regional-average daily radiative forcing (${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$) occurs in NEC, with values of ∼ 12 W m⁻², and the lowest regional average occurs in the Arctic, with values of ∼ 2.6 W m⁻². Regional-average radiative forcing for NA and EUA are both intermediate, with values of ∼ 3.1 and ∼ 3.5 W m⁻², respectively.

Table 1 Statistics for regionally averaged (5th and 95th quantiles) albedo reduction ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$) and daily radiative forcing (${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$; W m⁻²).

Figure 6 Spatial distributions of averaged (a) $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and (b) ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ and statistics for regionally averaged (c) $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and (d) ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ for the Northern Hemisphere ISCA in December–May during the period 2003–2018. The boxes denote the 25th and 75th quantiles, and the horizontal lines represent the 50th quantiles (medians); the averages are shown as red dots; the whiskers denote the 5th and 95th quantiles.

On a regional level, NEC $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ falls primarily within the range ∼ 0.077–0.14, and intra-regional variability is relatively small due to pervasive heavy pollution (Fig. S2). Compared to snow albedo reduction, the radiative forcing for NEC exhibits a slightly greater spatial variability due to latitude-dependent differences in the flux of surface solar radiances, ranging from ∼ 7.2 W m⁻² to ∼ 17 W m⁻². In NA, where the principal ISCAs are located in southern Canada, the western US, and northern US Great Plains, $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ tend to range between ∼ 0.014–0.046 and ∼ 1.3–7.0 W m⁻², respectively. In EUA, $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ fall largely within the respective ranges of ∼ 0.017–0.049 and ∼ 1.6–8.4 W m⁻². Central Asia and Mongolia exhibit relatively high values for $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ (> 0.04) and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ (> 2 W m⁻²), while this pattern likely reflects the influence of anthropogenic BC in addition to natural dust (Pu et al., 2017; Zhou et al., 2019) (Fig. S2a–b).

In the Arctic, $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ both present quite large intra-regional variabilities, from ∼ 0.0028 to ∼ 0.046 and ∼ 0.48 to 6.6 W m⁻². Greenland has the cleanest snow, with $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ of ∼ 0.011–0.023 and ∼ 0.40–3.3 W m⁻². In the Canadian Arctic, $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ are mainly in a range of ∼ 0.012–0.055 and ∼ 0.59–6.1 W m⁻². In addition, the relatively high values are found around the edge of ISCA to the west of the Canadian Arctic. The possible reason is that these areas are suffering from faster snow melting compared with the rest of the Canadian Arctic in spring, which is characterized by higher snow grain size (Fig. 4b). Hence, more LAPs are accumulated in the surface snow, resulting in higher snow albedo reduction. In the Russian Arctic, $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ values increase with altitude by ∼ 0.012–0.048 and ∼ 1.0–7.3 W m⁻². The snow albedo reduction in eastern Siberia is quite high and comparable with the values at midlatitudes. Moreover, benefiting from the higher solar radiances in eastern Siberia in spring (Fig. S2d) than those at midlatitudes in winter–spring (Figs. 4c and S2c), ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ in eastern Siberia is higher than parts of the midlatitudes. Even different from the findings in previous modeling studies (e.g., Flanner et al., 2007, 2009), the results seem to be comparable with the limited ground-based estimates (Fig. S3). The serious biomass burning in eastern Siberia in spring may be responsible for such high values (Warneke et al., 2010; Hegg et al., 2009). Overall, the Arctic spatial pattern of $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ in our study is consistent with the previous studies based on field experiments (Dang et al., 2017) and model simulation (Flanner et al., 2007). Nevertheless, we note that readers should be cautious about our reported high values in the Russian Arctic, and more field experiments are necessary for validating the results.

As mentioned above, the assumption of semi-infinite snowpack will trigger an overestimate for radiative forcing when snow depth is not thick enough. Figure 7 shows the spatial distribution of the ratio of retrieved radiative forcing using semi-infinite snow to radiative forcing using ERA-Interim snow depth. As can be seen, semi-infinite snowpack assumption will lead to an overestimate of up to ∼ 25 % in midlatitude areas, where snow depth is thin. In contrast, the influence of snow depth on radiative forcing is negligible in the Arctic, where snow is thick enough to become semi-infinite snowpack. These results demonstrated the important impact of snow depth on radiative forcing retrievals, which must be considered to reduce the overestimate for the following study.

Figure 7 The spatial distribution of the ratio of retrieved radiative forcing using semi-infinite snow to radiative forcing using ERA-Interim snow depth.

4.5 Accuracy discussion

In spite of the rigorous processes for radiative forcing retrieval, the uncertainty still existed. For example, light-absorbing particles in the atmosphere will reduce the accuracy of MODIS surface reflectance retrieval even though the atmospheric correction has been conducted. In addition, previous studies have pointed out a high scatter when converting NDSI to FSC using Eq. (5), which will induce bias in snow albedo retrieval (Rittger et al., 2013; Riggs et al., 2016; Bair et al., 2019). Furthermore, the method for snow grain size retrieval is only based on a single MODIS band at 1.24 µm, which could lead to higher uncertainties. Above all, all of these factors will result in a non-negligible uncertainty for radiative forcing retrieval, which needs to be further discussed.

To account for this issue, we consider that the accuracy of atmospheric correction is typically ±(0.005 + 0.05 × reflectance) under conditions under which aerosol optical depth (AOD) is less than 5.0, and solar zenith angle is less than 75^∘ according to the MODIS Surface Reflectance user's guide (Collection 6; https://modis.gsfc.nasa.gov/data/dataprod/mod09.php, last access: 1 October 2020). In addition, the bias for FSC calculation is assumed to be 10 % according to Riggs et al. (2016). The bias for snow grain size retrieval is assumed to be 30 % according to the studies of Pu et al. (2019) and Wang et al. (2017). Figure 8 shows the overall uncertainty in radiative forcing retrieval due to all these factors, while Fig. S6 shows the uncertainty caused by each factor. In general, the upper (lower) bound of the uncertainty falls within a range of 15 %–108 % (−106 % to −20 %), with atmospheric correction and FSC calculation contributing more to the uncertainty than snow grain size retrieval. The highest uncertainty occurs in the Arctic, while the lowest uncertainty occurs in NEC. Furthermore, the uncertainty shows a negative correlation with retrieved radiative forcing. The results indirectly demonstrated the reasonability of different correction factors performed in different regions. For example, the value of 1.6 used in NEC suggests that the correction approach works well for heavily polluted snow, while the value of 6.0 used in Greenland for relatively clean snow suggests that the method becomes not accurate enough.

Figure 8 The overall lower bound and upper bound of the uncertainty range of radiative forcing retrieval due to atmospheric correction, MODIS-derived snow grain size retrieval, and snow cover fraction calculation.

It is worth noting that the uncertainties from these factors could not fully explain the high correction factor in clean snow. The reason for why the ratio $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,daily}}^{\mathrm{LAPs}}/\mathrm{\Delta }{\mathit{\alpha }}_{\text{in situ,daily}}^{\mathrm{LAPs}}$ is larger than 1 is most likely the effect of snow surface roughness (Manninen et al., 2020) and vegetation (Pu et al., 2019). This effect is not accounted for in SNICAR and probably reduces the derived albedo from MODIS and therefore results in an overestimation of the albedo reduction attributed to LAPs. Moreover, there are other potential factors causing errors: (1) MODIS has variably spaced and discrete spectral bands and thus cannot provide a continuous spectral measurement of reflectance, resulting in a non-negligible uncertainty in retrieving the radiative forcing by LAPs in snow (Painter et al., 2012a); (2) we use the retrieved radiative forcing in a pixel size of 0.05^∘ × 0.05^∘ to compare with the in situ radiative forcing calculated from the measured BC_equiv concentration with a sample site located somewhere within the pixel; however, such a comparison may not be representative at some sites due to the inhomogeneous spatial distribution of LAP contents, which will influence radiative forcing retrieval; (3) in situ measurements also have uncertainties, which may cause a high bias for snow albedo reduction in clean snow. For example, a 10 % bias for 50 ng g⁻¹ BC can result in an 8 % bias for snow albedo reduction.

4.6 Attribution to the spatial variability in snow albedo reduction and radiative forcing

Here, we address the attributions to the spatial variability in snow albedo reduction and radiative forcing. As discussed in Sect. 3.2.5, the spatial variability in snow albedo reduction and radiative forcing is largely dependent on LAP content, snow grain radius, snow depth, and the geographic factor. Figure 9 illustrates the fractional contributions of each factor within the study regions. For the Northern Hemisphere ISCA as a whole, LAPs (I_LAPs) are the greatest contributor (84.3 %) to snow albedo reduction, followed by SD (13.7 %); R_eff and G have only a minor influence (1.9 % and < 1 %, respectively) (Fig. 9a). This result confirms that the concentration of LAPs in the snowpack plays a fundamental role in spatial variability in snow albedo reduction.

Figure 9 Fractional contributions of LAPs, snow grain size (R_eff), geographic factor (G), and snow depth (SD) to the spatial variations in (a) snow albedo reduction and (b) daily radiative forcing.

LAPs also constitute the dominant contributors to snow albedo reduction on a regional scale, accounting for 96.0 % of the Arctic signal, 56.7 % in EUA, and 49.9 % in NA, and are the second-largest contributor in NEC (40.3 %). The contribution of SD is greatest in NEC (56.3 %), with slightly lower values in EUA (40.3 %) and NA (48.8 %), reflecting the significant spatial variability in SD across these regions. In the Arctic, the snowpack is sufficiently thick to be considered a homogeneous, semi-infinite snowpack, and thus the contribution of SD is negligible. In contrast, R_eff makes only minor contributions in NEC (3.3 %), NA (1.3 %), EUA (2.8 %), and the Arctic (1.4 %). Finally, G makes the smallest contribution to snow albedo reduction (< 1 %), both on regional and global scales.

On a hemispheric scale, the greatest contributors to radiative forcing are LAP content (70.0 %) and G (22.3 %), followed by SD (7.6 %). As with snow albedo reduction, R_eff plays only a minor role. The influence of G on spatial variability in radiative forcing is attributed to the high degree of variability in latitude-dependent solar radiative fluxes among ISCA. On a regional scale, the respective contributions of LAP content, G, and SD are also comparable among the four study areas, accounting for 34.1 %, 11.1 %, and 52.0 % of radiative forcing in NEC; 39.2 %, 13.9 %, and 46.4 % in NA; and 48.0 %, 19.3 %, and 31.6 % in EUA. The Arctic radiative forcing is dominated by LAPs (85.6 %) and G (12.7 %).

In summary, LAPs play a dominant role in the spatial variability in snow albedo reduction and radiative forcing. Our results also highlight the significant contribution of SD to snow albedo reduction and of G to radiative forcing.

4.7 Comparisons with model simulations

To investigate the global distribution and variance of RFLS, previous studies have tended to rely on Earth system models with minimal cross-checking from in situ measurements or remote-sensing observations (Qian et al., 2015; Skiles et al., 2018). In this study, we compared MODIS retrievals with CESM2 to improve our understanding of the magnitude of RFLS on a global scale.

Employing snow BC concentrations from CESM2, we also calculated December–May daily radiative forcing (RF_CESM2) for the Northern Hemisphere ISCA during the period 2003–2014 (Fig. 10a). Statistics are presented in Fig. S7. Briefly, RF_CESM2 exhibits strong spatial inhomogeneity, with values ranging from 0.20 to 5.6 W m⁻². The highest regional average in RF_CESM2 occurs in NEC (≥10 W m⁻²) and the lowest in the Arctic (≤0.5 W m⁻²), consistent with ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$.

Figure 10 (a) Spatial distributions of daily average radiative forcing (RF_CESM2) based on the CESM2 soot content of snow in December–May for the period 2003–2014. (b) Scatterplot of ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ versus RF_CESM2.

Figure 10b depicts the comparison of ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ and RF_CESM2. In NEC, RF_CESM2 (15 W m⁻²) compares well with ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ (12 W m⁻²), with a significant correlation at the 99 % confidence level. For EUA, RF_CESM2 (3.8 W m⁻²) is similar to ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ (3.5 W m⁻²). For NA, RF_CESM2 (1.2 W m⁻²) is lower than ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ (3.1 W m⁻²), and the spatial correlation between them is poor. In the Arctic, RF_CESM2 is correlated with ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ at the 99 % confidence level. However, RF_CESM2 (1.7 W m⁻²) is lower than ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ (2.6 W m⁻²) by a factor of 1.5.

Overall, the RFLS derived from our MODIS retrievals and modeling-based estimates exhibits the same magnitude over the Northern Hemisphere. In NEC, the MODIS-derived and model-derived estimates show good general agreement, indicating the satisfactory performance of CESM2 in this heavily polluted region. In EUA, average radiative forcing values are comparable, but the spatial correlation is relatively poor, while MODIS retrievals for the Arctic are significantly higher than those simulations.

5 Discussion

In recent decades, there has been increasing scientific interest in snow LAPs due to their role in the climate system, and numerous studies have attempted to evaluate RFLS. In addition to making global-scale comparisons between our MODIS retrievals and model-based estimates, this study collects a comprehensive set of radiative forcing estimates based on local-scale observations and remote sensing to make quantitative regional- and global-scale comparisons and synthetically evaluate the magnitude of RFLS (Table 2). This approach also affords the opportunity to examine the MODIS retrievals used in our study.

Table 2 Comparisons of radiative forcing due to LAPs in snow (this study) with observed and model-simulated values from previous studies.

Dang et al. (2017) reported RFLS values of 7–18, 0.6–1.9, and 0.1–0.8 W m⁻² for northern China, North America, and the Arctic, respectively, which only focused on the period of January–March and therefore are smaller than our retrievals. In NA, Sterle et al. (2013) estimated a daily averaged RFLS of ∼ 2.5–40 W m⁻² for the eastern Sierra Nevada in February–May 2009, while Miller et al. (2016) reported a daily RFLS of ∼ 35–86 (37–100) W m⁻² based on in situ (remote-sensing) measurements in the San Juan Mountains in May 2010. Both values are higher than our estimate (∼ 3.1 W m⁻²), potentially due to the significant dust deposition in those areas.

We also collected the daily average RFLS simulated by regional and/or global climate models (Table 2). For NEC, Zhao et al. (2014) and Qian et al. (2014) reported values of 10 W m⁻² in January–February and 5–10 W m⁻² in April, respectively. In NA, Qian et al. (2009) provided an estimate of 3–7 W m⁻² for the central Rockies and southern Alberta in March, while Oaida et al. (2015) reported an average RFLS of 16 W m⁻² over the western US in spring. Finally, Qian et al. (2014) and Qi et al. (2017) estimated RFLS values of < 0.3 and 0.024–0.39 W m⁻² for the Arctic in April, respectively. We consider our retrievals for NEC to be comparable with these regional model simulations, despite some disparity. However, we note that our result is significantly lower than those of previous studies in NA but higher in the Arctic.

On a global scale, Hansen and Nazarenko (2004) reported an RFLS of 0.3 W m⁻², while Flanner et al. (2007) showed an RFLS of ∼ 0.05 W m⁻². For the Northern Hemisphere as a whole, Bond et al. (2013) estimated a climate forcing of 0.13 W m⁻². Each of these previous values is significantly lower than our retrieval (∼ 2.9 W m⁻²). However, those studies included all areas regardless of snow cover throughout the whole year, while our results are only for Northern Hemisphere ISCA from December to May.

Overall, we consider our MODIS-based retrievals to be physically realistic on both regional and global scales, although we note a number of differences between our results and those generated by different methods. On the other hand, while in situ measurements are the most precise, their spatial coverage is restricted by logistical limitations and the extreme environments involved. Conversely, models can provide broad perspectives of climatic impacts yet are typically undermined by large uncertainty. Therefore, we argue that remote sensing provides a powerful technique with high spatial and temporal resolutions that can bridge the gap between in situ measurements and climate models and reduce the uncertainties associated with the latter. Further retrieval of remote-sensing data, including the use of multiple satellites and sensors, is therefore warranted to exploit this opportunity fully. We also indicate the fact that in parts of central EUA and the Russian Arctic, however, studies are barely performed but desired. Finally, we note that in situ observations remain limited, and more field campaigns are needed to constrain remote-sensing retrievals and modeling simulations.

6 Conclusions

We presented a global-scale evaluation of the daily radiative forcing of LAPs in the Northern Hemisphere snowpack (RFLS), estimated from remote-sensing data. The satellite-retrieved RFLS also has implications for expanding the value of limited in situ measurements, which can provide valuable information for climate models and help optimize model simulations.

Based on the corrected snow albedo reduction ($\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$), we calculated daily average RFLS (${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$) during December–May for the period 2003–2018. For the identified snow-covered area over the Northern Hemisphere as a whole, average $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ is ∼ 0.021, and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ is ∼ 2.9 W m⁻². We also observed distinct spatial variability in snow albedo reduction and RFLS. The highest regional-average $\mathrm{\Delta }{\mathit{\alpha }}_{\text{MODIS,corrected}}^{\mathrm{LAPs}}$ (∼ 0.11) and ${\mathrm{RF}}_{\text{MODIS,daily}}^{\mathrm{LAPs}}$ (∼ 12 W m⁻²) occur in northeastern China, while the lowest regional averages of ∼ 0.016 and ∼ 2.6 W m⁻², respectively, are observed in the Arctic. Moreover, we indicated that the semi-infinite assumption could overestimate up to ∼ 25 % of RFLS, especially for thin and patchy snow, such as midlatitudes in Eurasia and NA. In addition, if the ground-based corrections are not considered, the total uncertainty in RFLS retrievals is in the range of 15 %–108 % (lower bound of −106 % to −20 %) due to atmospheric correction, snow cover fraction calculation, and snow grain size retrieval.

Following this assessment, we made quantitative attributions of the spatial variability in snow albedo reduction and radiative forcing. Our results indicate that the LAP content is the largest contributor (84.3 %) to spatial variance in snow albedo reduction, followed by snow depth (13.7 %), whereas snow grain size (1.9 %) and the geographic factor G (< 1 %) are only minor contributors on a northern-hemispheric scale. LAP content and G account for 70.0 % and 22.3 % of the spatial variability in radiative forcing, respectively, followed by SD (7.6 %) over the Northern Hemisphere.

Retrieved RFLS values are compared spatially with the model-derived estimates of the CESM2. Our results indicate that MODIS retrievals show the same magnitude with modeled estimates for the Northern Hemisphere. However, although the CESM2 performs well in NEC, there remain large uncertainties in the Arctic. To evaluate and examine the MODIS retrievals synthetically, we then compared the retrieved RFLS to previously published estimates, including local-scale observations, remote-sensing retrievals, and regional- and global-scale model simulations. The results of this evaluation suggest that MODIS retrievals are generally realistic, despite a number of important differences among the various methods.

Finally, we urge the community to expand the ground-based measurements of the global snowpack, particularly in those regions currently lacking in situ observations. Such development would help further constrain and improve satellite-based retrievals in the future. We propose that climate models validated by these refined remote-sensing retrievals should be able to capture the RFLS more accurately, thereby providing more reliable estimates of the future impacts of global climate change.