The remote sensing of radiative forcing by light-absorbing particles (LAPs) in seasonal snow over northeastern China

Light-absorbing particles (LAPs) deposited on snow can decrease snow albedo and affect climate through snow-albedo radiative forcing. In this study, we use MODIS observations combined with a snow-albedo model (SNICAR – Snow, Ice, and Aerosol Radiative) and a radiative transfer model (SBDART – Santa Barbara DISORT Atmospheric Radiative Transfer) to retrieve the instantaneous spectrally integrated radiative forcing at the surface by LAPs in snow (${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$) under clear-sky conditions at the time of MODIS Aqua overpass across northeastern China (NEC) in January–February from 2003 to 2017. ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ presents distinct spatial variability, with the minimum (22.3 W m⁻²) in western NEC and the maximum (64.6 W m⁻²) near industrial areas in central NEC. The regional mean ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is $\sim \mathrm{45.1}\pm \mathrm{6.8}$ W m⁻² in NEC. The positive (negative) uncertainties of retrieved ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ due to atmospheric correction range from 14 % to 57 % (−14 % to −47 %), and the uncertainty value basically decreases with the increased ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. We attribute the variations of radiative forcing based on remote sensing and find that the spatial variance of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in NEC is 74.6 % due to LAPs and 21.2 % and 4.2 % due to snow grain size and solar zenith angle. Furthermore, based on multiple linear regression, the BC dry and wet deposition and snowfall could explain 84 % of the spatial variance of LAP contents, which confirms the reasonability of the spatial patterns of retrieved ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in NEC. We validate ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ using in situ radiative forcing estimates. We find that the biases in ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ are negatively correlated with LAP concentrations and range from ∼5 % to ∼350 % in NEC.

1 Introduction

Pure snow is the most strongly reflective natural substance at the surface of the Earth, and seasonal snow covers more than 30 % of the Earth's land area (Painter et al., 1998). Therefore, snow cover has an important impact on the radiation balance of the Earth (Cohen and Rind, 1991). Light-absorbing particles (LAPs), such as black carbon (BC), organic carbon (OC), and mineral dust, deposited on snow can effectively reduce snow albedo (Hadley and Kirchstetter, 2012; He et al., 2017, 2018; Li et al., 2016; Warren, 1982, 1984; Warren and Wiscombe, 1980) and enhance the absorption of solar radiation (Dang et al., 2017; Kaspari et al., 2014; Liou et al., 2011, 2014; Painter et al., 2012b). Warren and Wiscombe (1980) determined that 10 ng g⁻¹ of BC in old snow could reduce the snow albedo by nearly 1 % at 400 nm with a snow grain size of 1000 µm. Jacobson (2004) pointed out that the snow-albedo reduction caused by BC in snow and ice is 0.4 % on a global scale and 1 % in the Northern Hemisphere based on model simulations. LAPs in snow further contribute to alterations in snow morphology, accelerations in snowmelt, and reductions in snow cover (Flanner et al., 2007, 2009; Painter et al., 2013a; Xu et al., 2009). For example, Qian et al. (2009) found that simulated BC-induced snow-albedo perturbations lead to a significant decrease in snow water equivalent by 2–50 mm over the mountains during late winter to early spring in the western United States. Ming at al. (2015) pointed out that widespread albedo decrease and the induced melting of Himalayan snow and ice in the early 21st century are partly caused by LAP deposition, resulting in approximately 10.4 Gt yr⁻¹ of mass loss equivalent of the Hindu Kush, Karakoram, and Himalaya (HKH) glaciers.

Several studies pointed out that the radiative forcing effects of LAPs in snow on local hydrological cycles (Painter et al., 2010; Qian et al., 2009; Yasunari et al., 2010) and regional and global climate (Bond et al., 2013; Hansen and Nazarenko, 2004; He et al., 2014; Jacobson, 2002, 2004; McConnell et al., 2007; Ramanathan and Carmichael, 2008; Yasunari et al., 2015) are non-negligible based on model simulations. In the Northern Hemisphere, Hansen and Nazarenko (2004) illustrated that the radiative forcing of BC on snow and ice albedo is +0.3 W m⁻². In addition, the IPCC AR5 (2013) indicated that the impact of BC in snow and ice accounted for a global mean climate forcing of +0.04 W m⁻², but the confidence level is low. Bond et al. (2013) estimated climate forcing consisting of radiative forcing, rapid adjustments, and strong snow-albedo feedback due to BC-in-snow forcing and pointed out that the best valuation of climate forcing by BC in snow and sea ice is +0.13 W m⁻², with the 90 % uncertainty bounds ranging from +0.04 to +0.33 W m⁻². Nevertheless, recent studies reported that ample factors confuse model simulations of BC-in-snow-induced climate forcing, and the model-based estimate of the regional and global radiative forcing caused by BC in snow and ice is still a challenge (Hansen and Nazarenko, 2004; Bond et al., 2013; Pu et al., 2017).

Much of northeastern China (NEC) is covered by contiguous seasonal snow in the winter and early spring. Local pollutant emissions in these regions are some of the most intense in the world (Bond et al., 2004), leading to considerable amounts of LAPs deposited into snow via wet and dry deposition (Bond et al., 2013). Therefore, several field campaigns have been conducted to investigate the LAP concentrations in snow across NEC (Huang et al., 2011; X. Wang et al., 2014, 2015). Wang et al. (2013) conducted a large field campaign to measure LAPs in seasonal snow across northern China from January to February 2010. They found that BC is the dominant absorber compared with OC and dust in NEC, and BC concentrations in seasonal snow range from 40 to 4000 ng g⁻¹, which are much higher than those measured in the Arctic, North America, and Europe (Doherty et al., 2010, 2014; Peltoniemi et al., 2015). Recently, Wang et al. (2017) showed that LAPs can reduce the visible spectral albedo by ∼0.35 in NEC based on in situ measurements and model simulations, which indicated a significant impact of LAPs on snow-albedo reduction. Zhao et al. (2014) simulated radiative forcing by LAPs in snow over northern China using a coupled model, and they noted that the uncertainties of their results are non-negligible due to limited observations.

Remote sensing is considered to be a powerful tool for estimating snow physical properties (e.g., Nolin and Dozier, 1993, 2000). Snow spectral albedo is highly dependent on wavelength λ. The albedo of pure snow is extremely high at visible (VIS) wavelengths (∼0.99 at λ=500 nm) but drops to a very low level in the near infrared (NIR; λ>1000 nm), where the imaginary part of the complex refractive index for ice is orders of magnitude greater than that in the VIS wavelengths (Wiscombe and Warren, 1980). The NIR albedo is sensitive to snow grain size; as grain size increases, the photon paths through ice get longer so there is a greater absorption probability. The NIR albedo is also sensitive to solar zenith angle: at low sun a photon's first scattering event occurs closer to the surface so it is more likely to escape (Wiscombe and Warren, 1980; Warren, 2013). Previous studies have successfully retrieved snow grain size using satellite NIR albedo data and radiative transfer models (e.g., Nolin and Dozier, 2000). On the other hand, the VIS albedo of snow is insensitive to grain size and solar zenith angle, which means that natural aging-induced changes in snow grains have little effect on VIS snow albedo. However, the VIS snow albedo is instead sensitive to LAPs in a semi-infinite snowpack. When LAPs such as BC or dust are present, snow albedo decreases primarily in the VIS wavelengths (Ming et al., 2012; Wang et al., 2017). This albedo reduction results from the greater imaginary part of the complex refractive index for LAPs compared with that of the highly transparent ice, which leads to more light absorption (Warren and Wiscombe, 1980). Therefore, the snow spectral albedo derived from satellite remote sensing in the VIS wavelengths can be used to estimate the impact of LAPs on snow albedo, which furthermore provides valuable information for modeling simulations to reduce relative uncertainties. To estimate the influence of mineral dust on snow albedo in the European Alps, Di Mauro et al. (2015) defined a new spectral index, the snow darkening index, based on in situ measured snow spectral reflectance and Landsat 8 Operational Land Imager (OLI) data; they found that the snow darkening index could effectively track the content of mineral dust in snow. In addition, Di Mauro et al. (2017) characterized the impact of LAPs on ice and snow albedo of the Morteratsch Glacier, a large-valley glacier in the Swiss Alps, using satellite (EO-1 Hyperion) hyperspectral data. The results showed that the spatial distribution of both narrowband and broadband indices retrieved from Hyperion was highly correlated with ice and snow impurities. In the Arctic, Dumont et al. (2014) developed an impurity index based on satellite observations (MODIS C5 surface reflectance) to analyze the snow darkening caused by the increased contents of LAPs in snow in Greenland. Nevertheless, Polashenski et al. (2015) pointed out that the apparent snow albedo declines in Greenland observed from MODIS C5 surface reflectance (Dumont et al., 2014) have a significant contribution from uncorrected Terra sensor degradation. In this study, in order to prevent interference from sensor degradation, we used the latest version (version 6, C6) of MODIS data from the Aqua sensor, which does not suffer from the influence of sensor degradation (Polashenski et al., 2015). Even though these studies have confirmed the ability of remote sensing to assess the role of LAPs in snow-albedo reduction, they did not quantitatively estimate the radiative forcing due to LAPs in snow, which is extremely important for inferring the impact of LAPs on regional and global climate. Recently, Ming et al. (2012) estimated the radiative forcing in Himalayan glaciers based on the differences between the simulated pristine albedo and the satellite observation albedo, which could be partly attributed to BC and dust. The results illustrated that the current surface radiation absorption could lead to significant melting in Himalayan glaciers, which could cause most of them to be in danger of rapid mass loss. Furthermore, Painter et al. (2012a) successfully used the MODIS Dust Radiative Forcing in Snow (MODDRFS) model to retrieve surface radiative forcing by LAPs in snow cover from Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance data. They found that instantaneous at-surface radiative forcing can be beyond 250 W m⁻² in the Hindu Kush–Himalaya area and falls in a range of 30–250 W m⁻² in the upper Colorado River Basin. Painter et al. (2013b) also provided and validated an algorithm suite to quantitatively retrieve radiative forcing by LAPs in snow from Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data in the Senator Beck Basin Study Area (SBBSA), SW Colorado, USA. The lowest radiative forcing was found on the high north-facing slopes, while the highest was on southeast-facing slopes at the lowest elevations. Seidel et al. (2016) analyzed the spatial and temporal distribution of radiative forcing by LAPs in snow in the Sierra Nevada and Rocky Mountains from imaging spectroscopy. Their results presented an increased radiative forcing from 20 up to 200 W m⁻² in the melting period. Warren (2013) also indicated that attempts to use satellite remote sensing to estimate the radiative forcing by LAPs in polluted regions are likely feasible. However, to date, no studies have quantitatively investigated the contributions of each factor to the variations of radiative forcing by LAPs in snow based on remote sensing. Moreover, the radiative forcing by LAPs in snow across NEC is far less studied by using satellite remote sensing, even though the LAP contents in these regions are much higher compared with those in the Arctic, Europe, and the USA (Dang et al., 2017).

Estimating the radiative forcing by LAPs in snow by using surface measurements is more precise than using remote sensing or model simulation. However, surface measurements of snow albedo and LAP content in snow are very limited on regional or global scales. Until now, the observational sample sites (<50) have been sparse and only for two individual measurements in 2010 and 2014 over a wide NEC area of ∼1.5 million km² (X. Wang et al., 2013, 2017; Z. W. Wang et al., 2014; Ren et al., 2017). The very limited number of measurement sites led to poor spatial–temporal distribution of the estimated radiative forcing in NEC (Dang et al., 2017). On the other hand, remote sensing technology has the advantage of high spatial–temporal resolution and has been successfully used to retrieve the radiative forcing by in-snow light-absorbing particles in high-snow-cover areas (Painter et al., 2012a). In addition, a previous study indicated that the uncertainty in estimating radiative forcing using model simulations is very high due to limited measurement data (Zhao et al., 2014); however, this could possibly be improved by combining remote-sensing-retrieved results. Hence, estimating the radiative forcing by LAPs in snow by using satellite remote sensing seems to be necessary.

In this study, we attempt to retrieve the radiative forcing by LAPs in snow across NEC using MODIS datasets combined with the Snow, Ice, and Aerosol Radiation (SNICAR) model (Flanner et al., 2007, 2009) and the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) model (Ricchiazzi et al., 1998); we estimate the uncertainties of radiative forcing from atmospheric correction and qualify the fractional contribution of each factor to the spatial variance of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. Then, we investigate the reasonability of the spatial patterns of retrieved radiative forcing in NEC based on BC deposition and snowfall data. Finally, we quantitatively estimate the biases of MODIS-retrieved radiative forcing using in situ radiative forcing estimates, which are based on snowfield measurements.

2 Datasets

2.1 Remote sensing datasets

The latest version (Collection 6) of MODIS surface reflectance data (MYD09GA), MODIS snow cover data (MYD10A1), and MODIS aerosol optical depth (AOD) data (MYD04) is used in this study from 2003 to 2017, which covers the months of January through February (https://modis.gsfc.nasa.gov/, last access: 20 March 2018). The MOD09 product is divided into seven bands (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; and band 7, 2105–2155 nm) and has a spatial resolution of 500 m (Vermote, 2015). The MOD09 surface reflectance is an estimate of the surface spectral reflectance for each band, which corrects for the effects of atmospheric gases and aerosols. The performance of the atmospheric correction algorithm suffers from the influence of view and solar zenith angles and aerosol optical thickness; the accuracy of the algorithm is also affected by the wavelengths of different bands. More details about the data products and a band quality description of MOD09GA can be found in the MODIS Surface Reflectance User Guide (https://modis.gsfc.nasa.gov/data/dataprod/mod09.php, last access: 20 March 2018). MODIS satellite data have been widely accepted in retrievals of snow cover and its physical properties (e.g., Scambos et al., 2007; Rittger et al., 2013). In addition, MODIS has three visible bands (VIS) and the radiometric range in the VIS over snow surface has no saturation phenomenon, which provides the ability to detect changes in reflectance in the VIS caused by LAPs in snow (Painter et al., 2012a).

2.2 Surface measurement datasets

Wang et al. (2017) conducted a snow survey across NEC in January 2014. They measured AOD using a Microtops II sun photometer. The Microtops II sun photometer is a portable instrument and measures solar radiance in five spectral wave bands (340, 440, 675, 870, and 936 nm) from which it automatically derives aerosol optical depth (AOD). When the Microtops II sun photometer is well cleaned and well calibrated, its AOD retrievals can be comparable with those of CIMEL sun photometers used in the AERONET network, with uncertainties ranging from 0.01 to 0.02 (Ichoku et al., 2002). The snow albedo and surface solar irradiance were measured using an Analytical Spectral Devices (ASD) spectroradiometer. The ASD spectroradiometer has 3 nm spectral resolution on the visible–near-infrared detector (350–1050 nm, silicon photodiode array) and 10–12 nm resolution on the shortwave infrared detectors (900–2500 nm, InGaAs). Measurements are made by standing “down-sun” of the receptor, taking consecutive scans of downwelling and upwelling radiation. Wuttke et al. (2006) indicated that the ASD spectroradiometer is considered to be the most mobile, capable, and rapid for measuring spectral albedo during short time periods, especially in very cold regions. The cosine error is less than 5 % for solar zenith angles below 85^∘ at a wavelength of 320 nm. We use these datasets to validate the snow grain size retrievals and the simulated surface solar irradiance values.

Snow samples were collected at 46 sites in January and February 2010 across northern China (Wang et al., 2013) and at 13 sites in January 2014 across northeastern China (Wang et al., 2017). A detailed description of the procedures of snow collection and filtration has been presented by previous studies (Doherty et al., 2010, 2014; Wang et al., 2013). Briefly, in order to ensure that the collected snow samples are regionally representative and minimize the influence of local emission sources, sample locations are usually chosen at least 1 km upwind of the approach roads and railways and more than 50 km from cities and towns. In addition, efforts are made to collect samples in open areas in order to prevent the contaminations from the detritus of bushes and trees. Generally, snow samples are collected within a vertical resolution varied from ∼2 to 10 cm and usually at typically vertical intervals of 5 cm from the top to the bottom throughout the snowpack depth at each site. In the case of a visibly distinct layering, such as newly fallen snow at the surface layer or a melt layer, the snow at that layer is gathered individually. Right and left snow samples of two side-by-side vertical profiles are collected within each layer to make a comparison and average the snow sample pairs. All snow samples are maintained frozen to prevent the melting snow from influencing the LAP content. Usually every 3 to 4 d, snow samples are filtered at temporary laboratories set up in hotels. Simply, snow samples are melted and filtered through Nuclepore filters of 0.4 µm pore size. The samples “before” and “after” filtration are gathered and refrozen for the following chemical analysis, and the filters are used for optical analysis.

An integrating sphere–integrating sandwich spectrophotometer (ISSW) is applied to analyze the filters and quantify the spectral light absorption by LAPs in snow. ISSW was first described by Grenfell et al. (2011), modified by Wang et al. (2013) and Doherty et al. (2014), and has been used by some previous studies (Dang and Hegg, 2014; Pu et al., 2017; Zhou et al., 2017). Schwarz et al. (2012) have confirmed the performance of ISSW in quantifying LAP concentrations in snow by comparing with a single-particle soot photometer (SP2), although both SP2 and ISSW may suffer from non-negligible uncertainties. Briefly, ISSW produces a diffuse radiation field when white light illumination is transmitted into an integrating sphere; then the diffuse radiation passes through the filter from below and is measured by a spectrometer. By measuring a sample filter and a blank filter, ISSW acquires the light attenuation spectrum due to the loadings on the sample filter (Grenfell et al., 2011). Because the design is such that the measured filter is sandwiched between two integrating spheres, the light attenuation is nominally due to the absorption of LAPs on the filter, and the influence of light scattering is negligible (Doherty et al., 2014). ISSW measures light attenuation from 400 to 700 nm by benefitting from an optimal signal-to-noise ratio and then extends the full spectrum to a range of 350 to 750 nm by extrapolation (Pu et al., 2017). Calibration is done by measuring a set of fullerene (synthetic BC; Alfa Aesar, Inc., Ward Hill, MA, USA) filters with a range of known loadings. Then, the light attenuation spectrum of the sample filter is transformed into an equivalent BC mass loading by comparing against the standard filters. With the loaded area on the filter and the volume of filtered snow water, equivalent BC mass is converted to equivalent BC concentration (BC_equiv). In this study, we will use BC_equiv for all LAPs to calculate the in situ radiative forcing.

2.3 BC deposition and emission data

BC deposition has important effects on the radiative forcing by LAPs in snow (Seidel et al., 2016). Higher BC deposition indicates that greater amounts of BC are deposited on snow, reducing the snow albedo. To our knowledge, there are no measurement data for the spatial distribution of BC deposition in NEC. Therefore, we collected reanalysis data on BC deposition from the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2) in January–February from 2003 to 2017 and modeling data on BC deposition from the Coupled Model Intercomparison Project Phase 6 (CMIP6, the latest CMIP phase), including CESM2, CESM2-WACCM, and CNRM-ESM2-1 historical experiments in January–February from 2003 to 2014 (Eyring et al., 2016). So far, only the above three models in CMIP6 have provided BC deposition data. In our study, we prefer to use MERRA-2 data because these data are the latest atmospheric reanalysis data of the modern satellite era produced by NASA's Global Modeling and Assimilation Office (GMAO). They assimilate aerosol observations and other observation types to provide a viable ongoing climate analysis. The provided observable parameters and aerosol diagnostics have wide potential applications ranging from air quality forecasting to aerosol–climate interactions (Bocquet et al., 2015; Randles et al., 2016, 2017). In addition, the period range of MERRA-2 BC deposition data satisfies our study period of 2003–2017, but the CMIP6 data are only updated to 2014. We note that the results and conclusions based on different BC deposition data are similar (see Sect. 4.3).

Local BC emission density can also indicate the LAP contents in snow. Among all available BC emission density data, we use the data from the research group at Peking University (http://inventory.pku.edu.cn/home.html, last access: 15 June 2018; R. Wang et al., 2014) after taking spatial and temporal resolution, data period, data quality, and other factors into account. The BC emission density data we used are for January–February from 2003 to 2014 because they are only updated to 2014.

2.4 Snowfall and snow parameter data

Seidel et al. (2016) pointed out that snowfall can affect the radiative forcing by LAPs in snow. A higher frequency of snowfall implies that greater amounts of fresh snow, which has smaller snow grains than aged snow, are present at the surface, increasing the snow albedo (Z. W. Wang et al., 2014). In this study, we collected four types of snowfall data in January–February from 2003 to 2017, including surface observational data from the China Meteorological Administration (126 observation stations), the ERA-Interim reanalysis (http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/, last access: 25 April 2018), the Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2), and the National Centers for Environmental Prediction (NCEP) Climate Prediction Center (CPC) (https://www.esrl.noaa.gov/psd/data/gridded/data.cpc.globalprecip.html, last access: 12 April 2019). Figure S1 in the Supplement shows the spatial distribution of the observational stations over northeastern China. We note that the observation stations are limited in our study areas. Compared with the observed snowfall data, we also assessed snowfall data from the ERA-Interim reanalysis, the MERRA-2 reanalysis, and CPC in NEC. We found that the ERA-Interim reanalysis data are more consistent with surface observations (Fig. S2). Therefore, we prefer to use ERA-Interim for snowfall data in this study. But as with BC deposition data, the results and conclusions based on different snowfall data are similar (see Sect. 4.3).

To briefly describe the snow cover condition in NEC in January–February, we collect multiple types of snow parameter data including snow cover data (MYD10CM and MYD10C2) from MODIS products (https://modis.gsfc.nasa.gov/data/dataprod/mod10.php, last access: 9 April 2019), snow depth data from the Canadian Meteorological Centre (CMC) (https://nsidc.org/data/NSIDC-0447/versions/1, last access: 9 April 2019), and snow water equivalent data (GlobSnow-2) from European Space Agency (ESA) Global Snow Monitoring for Climate Research (http://www.globsnow.info/, last access: 9 April 2019).

3 Methods

3.1 Models

3.1.1 SNICAR model

Flanner et al. (2007) have presented a comprehensive description for the Snow, Ice, and Aerosol Radiative (SNICAR) model, which is the most widely used multilayer snow-albedo model in the field of atmospheric sciences. Here, we just briefly give a summary of SNICAR. SNICAR simulates radiative transfer in the snowpack based on the theory of Wiscombe and Warren (1980) and the two-stream multilayer radiative approximation of Toon et al. (1989). The input optical parameters (mass extinction coefficient, single-scattering albedo, and asymmetry factors) of snow grains and LAPs are offline calculated using Mie theory. In addition, the types of surface spectral distribution (clear or cloudy sky) and incident radiation (direct or diffuse) can be chosen by users, and users must specify the solar zenith angle if the incident flux is direct. In general, users should input the parameters involving the type of surface spectral distribution and incident radiation, number of snow layers, snow thickness, density, snow grain radius, and the type and concentration of LAPs in each snow layer, as well as the albedo of the underlying ground. Following a previous study (Painter et al., 2012a), we assume one-layer semi-infinite snow to drive the SNICAR model in this study.

3.1.2 SBDART model

In this study, we use the Santa Barbara DISORT Atmospheric Radiative Transfer (SBDART) model (Ricchiazzi et al., 1998) to simulate the surface solar irradiance. SBDART is one of the most widely used models to calculate the radiative transfer at the Earth's surface and within the atmosphere in both clear and cloudy sky. SBDART is a combination of the DISORT (Discrete Ordinate Radiative Transfer) radiative transfer module (Stamnes et al., 1988), low-resolution atmospheric transmission models, and Mie theory. The radiative transfer equations for a plane-parallel, vertically inhomogeneous, non-isothermal atmosphere numerically integrated in SBDART are based on DISORT, and light scattering by water droplets and ice crystals is from Mie theory. SBDART already considers all important processes that affect the ultraviolet, visible, and infrared radiation fields. The key components of SBDART include standard atmospheric models, cloud models, extraterrestrial source spectra, gas absorption models, standard aerosol models, and surface models. SBDART is well suited for widespread use in atmospheric radiation and remote sensing studies. More details about the SBDART model can be found in the paper of Ricchiazzi et al. (1998).

3.2 Retrieval methods

In this study we use BC as a representative to describe the effect of LAPs on snow albedo. Figure 1a shows the spectral snow albedo from 300 to 1400 nm. Gray areas show the typical spectral solar irradiance at the time of MODIS Aqua overpass (local time of 13:30) in January–February in NEC; the yellow column bars represent MODIS band passes. We can see that when LAPs such as BC are deposited on snow, they can effectively reduce snow albedo in the visible bands, which contain about half of the total solar radiation. For a snowpack with a snow grain radius of 100–300 µm, 100 ng g⁻¹ of BC in snow (a typical BC concentration in snow of the remote clean areas in NEC) can reduce snow albedo by ∼0.05–0.08 at 500 nm; 1000 ng g⁻¹ of BC in snow (a typical BC concentration in snow of the polluted industrial areas in NEC) can reduce snow albedo by ∼0.12–0.2. On the other hand, the effects of BC decrease at longer wavelengths in the near infrared (NIR). Moreover, when wavelengths exceed 1150 nm, snow albedo is dominated by the snow optical effective radius (R_eff) and is independent of LAPs. As shown in Fig. 1b, snow-albedo reduction is not only dependent on LAPs in snow, but also on snow grain size and solar zenith angle (θ). Generally, the reduction in snow albedo caused by BC increases with BC concentration and R_eff, whereas it decreases with the solar zenith angle (θ). Based on these characteristics, we retrieve R_eff, the reduction in snow albedo, and the radiative forcing by LAPs in this section.

Figure 1 (a) The spectral albedo of snow with different R_eff values and BC contents simulated using SNICAR. The column bars represent MODIS bands, and the gray areas represent the typical solar irradiance in winter in NEC. (b) The reduction in the 300–1240 nm spectral-weighted integrated snow albedo as a function of BC for different R_eff values and solar zenith angles (θ) is simulated using SNICAR. (c) The variations in the impurity index (I_LAPs) with BC content are simulated using SNICAR.

3.2.1 Snow cover

Three methods have been widely used in mapping snow-covered area using MODIS data. In the first method, “binary” maps, pixels are classified as either “snow free” or “snow covered” (Hall et al., 1995). However, significant errors exist in such maps, as pixels with a resolution of 500 m are not always completely covered by snow. The second method, the MODSCAG retrieval algorithm, is a fractional snow algorithm that is based on spectral mixture analysis (Painter et al., 2009). However, it cannot be applied in NEC due to limited information on the spectral reflectances of the vegetation, soils, and rock in this region. Therefore, we use the third method, which is based on the reflectances in the visible and NIR bands and the normalized difference snow index (NDSI):

$$\begin{array}{}\text{(1)}& \mathrm{NDSI}={\displaystyle \frac{R_{\mathrm{band}\mathrm{4}}-R_{\mathrm{band}\mathrm{6}}}{R_{\mathrm{band}\mathrm{4}}+R_{\mathrm{band}\mathrm{6}}}},\end{array}$$

where R_band4 and R_band6 are the surface reflectances in bands 4 and 6. Following Negi and Kokhanovsky (2011), an area is determined to be snow covered if the NDSI and the reflectance in band 4 both exceed 0.6. We note that the following analyses are only done over the defined snow-covered areas and periods.

3.2.2 Retrieval of snow grain size

Many methods have been used to retrieve snow grain size (e.g., Lyapustin et al., 2009; Nolin and Dozier, 1993). However, in NEC, the efficacy of most of these methods is limited, as the reflectances in bands 1–4 are seriously affected by LAPs in polluted snow (Fig. 1a), and the reflectances in bands 6–7 are not sensitive to R_eff. Hence, R_eff is retrieved at a wavelength of 1240 nm (the central wavelength of band 5) using SNICAR (Wang et al., 2017).

We validate the retrieved R_eff values using in situ measurements. The mean absolute error (MAE) is 71 µm, which is slightly higher than that reported by Painter et al. (2009). Nevertheless, the results are still credible because this study investigates a larger spatial scale than the previous study.

3.2.3 Impurity index

To assess LAP contents in snow, we use the surface reflectances in bands 4–5 to derive an impurity index (I_LAPs):

$$\begin{array}{}\text{(2)}& I_{\mathrm{LAPs}}={\displaystyle \frac{\ln \left(R_{\mathrm{band}\mathrm{4}}\right)}{\ln \left(R_{\mathrm{band}\mathrm{5}}\right)}}.\end{array}$$

This quantity increases with the LAP content but is almost independent of R_eff and θ (Fig. 1c). Di Mauro et al. (2017) have successfully used I_LAPs to assess the variations of LAP contents in the snow of the Morteratsch Glacier in the Swiss Alps. In this study, we did not retrieve the concentrations of LAPs because such a retrieval is constrained by many unknown factors, such as size distribution, optical properties, and the mixing state of LAPs (He et al., 2017, 2018; Painter et al., 2013a; Pu et al., 2017). Therefore, the conversion from satellite spectra to ground concentrations of LAPs will cause significant errors.

3.2.4 Retrieval of radiative forcing by LAPs in snow

Instantaneous surface solar irradiance at the time of MODIS overpass in January–February is simulated using the SBDART model (Ricchiazzi et al., 1998) with MODIS AOD data as inputs. Wang et al. (2017) have validated the MODIS AOD data using in situ measurements in NEC. For the other inputs, the typical values for midlatitude winter provided by SBDART are used. As a result, the normalized mean bias (NMB) is less than 2 % (Fig. S3).

We estimate the instantaneous spectrally integrated radiative forcing at the surface by LAPs in snow (${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}})$ under clear-sky conditions at the time of MODIS Aqua overpass, which is a function of solar irradiance and the difference between the MODIS spectral reflectance and a simulated clean-snow ($R_{\mathit{\lambda }}^{\mathrm{clean}\text{-}\mathrm{snow}})$ reflectance (Miller et al., 2016). $R_{\mathit{\lambda }}^{\mathrm{clean}\text{-}\mathrm{snow}}$ is simulated using the SNICAR model based on the retrieved R_eff and MODIS-derived solar zenith angle (θ). On the other hand, for MODIS spectral reflectance, because MODIS provides only discrete reflectances, we simulate a continuous spectral reflectance by fitting SNICAR to the MODIS data and derive the fitting parameters by minimizing the RMSE (Painter et al., 2009):

$$\begin{array}{}\text{(3)}& \mathrm{RMSE}={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\sum _{\mathit{\lambda }=\mathrm{band}\mathrm{1}}^{\mathrm{band}\mathrm{4}}{\left(R_{\mathit{\lambda }}^{\mathrm{model}}-R_{\mathit{\lambda }}^{\mathrm{MODIS}}\right)}^{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where RMSE is the root mean square error, and $R_{\mathit{\lambda }}^{\mathrm{model}}$ and $R_{\mathit{\lambda }}^{\mathrm{MODIS}}$ represent the simulated and MODIS-derived reflectances at a wavelength λ. Thus, ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is expressed as follows:

$$\begin{array}{}\text{(4)}& {\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}=\sum _{\mathit{\lambda }=\mathrm{300}\mathrm{nm}}^{\mathrm{1240}\mathrm{nm}}{E}_{\mathit{\lambda }}\cdot D_{\mathit{\lambda }}\cdot \mathrm{\Delta }\mathit{\lambda },\end{array}$$

where E_λ is the solar irradiance at a wavelength λ simulated by the SBDART model; D_λ is the difference between the clean-snow ($R_{\mathit{\lambda }}^{\mathrm{clean}\text{-}\mathrm{snow}})$ and simulated reflectances ($R_{\mathit{\lambda }}^{\mathrm{model}})$ at a wavelength λ, and Δλ is 10 nm.

3.2.5 Uncertainties

The uncertainties in radiative forcing retrievals are primarily due to terrain, liquid snow water, snow patchiness, protrusion of vegetation, and atmospheric correction. The study areas are located on smooth plains, and the content of liquid snow water is limited in the study regions in January and February (Wang et al., 2013). Moreover, both experimental and theoretical evidence shows that the effect of liquid water in snow on snow reflectance is small in the shortwave part of the spectrum but obvious at wavelengths of 0.95 and 1.15 µm (O'Brien and Munis, 1975; O'Brien and Koh, 1981; Wiscombe and Warren, 1980), which are not included in the MODIS bands used in our study. As a result, the effect of liquid water in snow on the calculations of snow grain size, I_LAPs, and radiative forcing is limited. Therefore, the effects of terrain and liquid snow water on MODIS retrievals could be negligible.

In our study, snow-covered area is determined if the NDSI and the reflectance in band 4 both exceed 0.6, which means that fractional snow cover (FSC) is larger than 0.86 according to the FSC equation ($\mathrm{FSC}=-\mathrm{0.01}+\mathrm{1.45}\cdot \mathrm{NDSI}$) from the MODIS Snow Products Collection 6 User Guide (http://nsidc.org/data/MYD10A1, last access: 20 March 2018). In January and February, snow depth is high and reaches its maximum depth in NEC; snow patchiness in high-snow-cover areas is mostly due to the protrusion of vegetation according to the observations of field campaigns (X. Wang et al., 2013, 2014). Thus, the MODIS-derived surface reflectance ($R_{\mathit{\lambda }}^{\mathrm{MODIS}}$) in a pixel of our study areas is not snow reflectance but a mixture of snow and vegetation reflectance. Therefore, we need to correct the errors of snow reflectance caused by the protrusion of vegetation. According to Painter et al. (2009), $R_{\mathit{\lambda }}^{\mathrm{MODIS}}$ could be expressed as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{R}_{\mathit{\lambda }}^{\mathrm{MODIS}}& ={\displaystyle \frac{E_{\mathit{\lambda }}\cdot \mathrm{FSC}\cdot R_{\mathrm{snow}}^{\mathit{\lambda }}+E_{\mathit{\lambda }}\cdot (\mathrm{1}-\mathrm{FSC})\cdot R_{\mathrm{vegetation}}^{\mathit{\lambda }}}{E_{\mathit{\lambda }}}}\\ & =\mathrm{FSC}\cdot R_{\mathrm{snow}}^{\mathit{\lambda }}+(\mathrm{1}-\mathrm{FSC})\cdot R_{\mathrm{vegetation}}^{\mathit{\lambda }},\end{array}\end{array}$$

where $R_{\mathit{\lambda }}^{\mathrm{MODIS}}$ is MODIS-derived surface reflectance at a wavelength λ, and E_λ is solar irradiance at a wavelength λ. FSC is the fractional snow cover, which could be derived according to the FSC equation. $R_{\mathrm{snow}}^{\mathit{\lambda }}$ and $R_{\mathrm{vegetation}}^{\mathit{\lambda }}$ represent snow and vegetation reflectance, respectively, at a wavelength λ. $R_{\mathrm{vegetation}}^{\mathit{\lambda }}$ is from the study of Siegmund and Menz (2005). Then $R_{\mathrm{snow}}^{\mathit{\lambda }}$ could be expressed as

$$\begin{array}{}\text{(6)}& R_{\mathrm{snow}}^{\mathit{\lambda }}={\displaystyle \frac{(R_{\mathit{\lambda }}^{\mathrm{MODIS}}-\left(\mathrm{1}-\mathrm{FSC}\right)\cdot R_{\mathrm{vegetation}}^{\mathit{\lambda }})}{\mathrm{FSC}}}.\end{array}$$

Finally, the accuracy of MODIS surface reflectance (MYD09GA) due to atmospheric correction is typically calculated based on the MODIS Surface Reflectance User Guide (Collection 6; https://modis.gsfc.nasa.gov/data/dataprod/mod09.php, last access: 20 March 2018) as follows:

$$\pm (\mathrm{0.005}+\mathrm{0.05}\cdot \mathrm{reflectance}),$$

which is suitable under conditions that AOD is less than 5.0 and θ is less than 75^∘. Therefore, we also estimate the uncertainty of MODIS retrievals from atmospheric correction. Briefly, the MODIS-derived snow reflectance ($R_{\mathrm{snow},\mathrm{uncertainty}}^{\mathit{\lambda }}$), which takes into account the accuracy of the atmospheric correction, is expressed as

$$\begin{array}{}\text{(7)}& R_{\mathrm{snow},\mathrm{uncertainty}}^{\mathit{\lambda }}=R_{\mathrm{snow}}^{\mathit{\lambda }}\pm (\mathrm{0.005}+\mathrm{0.05}\cdot R_{\mathrm{snow}}^{\mathit{\lambda }}),\end{array}$$

and the fractional uncertainty of MODIS-retrieved snow grain size (${\mathrm{FU}}_{R_{\mathrm{eff}}}$) could be expressed as

$$\begin{array}{}\text{(8)}& {\mathrm{FU}}_{R_{\mathrm{eff}}}={\displaystyle \frac{R_{\mathrm{eff},\mathrm{uncertainty}}-R_{\mathrm{eff}}}{R_{\mathrm{eff}}}},\end{array}$$

where R_{eff, uncertainty} is the SNICAR-simulated snow grain size using the snow reflectance of $R_{\mathrm{snow},\mathrm{uncertainty}}^{\mathrm{1240}}$. Similar to snow grain size, the fractional uncertainty of I_LAPs (${\mathrm{FU}}_{I_{\mathrm{LAPs}}}$) and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ (FU_RF) is

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathrm{FU}}_{I_{\mathrm{LAPs}}}={\displaystyle \frac{I_{\mathrm{LAPs},\mathrm{uncertainty}}-I_{\mathrm{LAPs}}}{I_{\mathrm{LAPs}}}},\text{(10)}& {\displaystyle }{\mathrm{FU}}_{\mathrm{RF}}={\displaystyle \frac{{\mathrm{RF}}_{\mathrm{MODIS},\mathrm{uncertainty}}^{\mathrm{LAPs}}-{\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}}{{\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}}}.\end{array}$$

We note that the positive and negative uncertainty is asymmetric due to the nonlinearity of the SNICAR model.

3.2.6 Attribution of the spatial variance of radiative forcing by LAPs in snow

As discussed above, ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is dependent on I_LAPs (the indicator of LAPs), R_eff, and θ and could be expressed as

$$\begin{array}{}\text{(11)}& {\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}=f(I_{\mathrm{LAPs}},R_{\mathrm{eff}},\mathit{\theta }).\end{array}$$

As a result, the spatial patterns of I_LAPs, R_eff, and θ determine the spatial pattern of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. First, we keep R_eff and θ spatially constant with values of the spatial averages (${\stackrel{\mathrm{‾}}{R}}_{\mathrm{eff}}$ and $\stackrel{\mathrm{‾}}{\mathit{\theta }})$. The radiative forcing is only dependent on the distribution of I_LAPs:

$$\begin{array}{}\text{(12)}& {\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)=f(I_{\mathrm{LAPs}},\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}},\stackrel{\mathrm{‾}}{\mathit{\theta }}).\end{array}$$

Similarly, we could obtain two additional equations:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)=f(\stackrel{\mathrm{‾}}{I_{\mathrm{LAPs}}},R_{\mathrm{eff}},\stackrel{\mathrm{‾}}{\mathit{\theta }}),\text{(14)}& {\displaystyle }{\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(\mathit{\theta }\right)=f(\stackrel{\mathrm{‾}}{I_{\mathrm{LAPs}}},\stackrel{\mathrm{‾}}{R_{\mathrm{eff}}},\mathit{\theta }).\end{array}$$

Then ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is fitted with ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)$, ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)$, and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(\mathit{\theta }\right)$ using multiple linear regression, and the fitted radiative forcing (${\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}$) is expressed as

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}& =a+b\cdot {\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)+c\cdot {\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)\\ & +d\cdot {\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(\mathit{\theta }\right),\end{array}\end{array}$$

where a, b, c, and d are regression coefficients. In our study, we find that ${\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}$ could explained 99.9 % of the variance of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ (Fig. S4). Therefore, we can attribute the variance of ${\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}$ instead of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ to estimate the fractional contribution of I_LAPs, R_eff, and θ to radiative forcing. Equation (15) can be written as

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}& -\stackrel{\mathrm{‾}}{{\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}}=\\ & b\cdot \left({\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)-\stackrel{\mathrm{‾}}{{\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}}\left(I_{\mathrm{LAPs}}\right)\right)\\ & +c\cdot \left({\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)-\stackrel{\mathrm{‾}}{{\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}}\left(R_{\mathrm{eff}}\right)\right)\\ & +d\cdot \left({\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(\mathit{\theta }\right)-\stackrel{\mathrm{‾}}{{\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}\left(\mathit{\theta }\right)}\right),\end{array}\end{array}$$

where ${\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}-\stackrel{\mathrm{‾}}{{\mathrm{RF}}_{\mathrm{Fit}}^{\mathrm{LAPs}}}$ is the radiative forcing anomaly (${\mathrm{RF}}_{\mathrm{Fit},\mathrm{anomaly}}^{\mathrm{LAPs}}$). Then, Eq. (16) can be written as

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{\mathrm{RF}}_{\mathrm{Fit},\mathrm{anomaly}}^{\mathrm{LAPs}}& =b\cdot {\mathrm{RF}}_{\mathrm{MODIS},\mathrm{anomaly}}^{\mathrm{LAPs}}\left(I_{\mathrm{LAPs}}\right)\\ & +c\cdot {\mathrm{RF}}_{\mathrm{MODIS},\mathrm{anomaly}}^{\mathrm{LAPs}}\left(R_{\mathrm{eff}}\right)\\ & +d\cdot {\mathrm{RF}}_{\mathrm{MODIS},\mathrm{anomaly}}^{\mathrm{LAPs}}\left(\mathit{\theta }\right).\end{array}\end{array}$$

According to Huang et al. (2016) and Huang and Yi (1991), the fractional contribution of I_LAPs to the variance of radiative forcing (${\mathrm{FC}}_{I_{\mathrm{LAPs}}}$) can be expressed as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& {\mathrm{FC}}_{I_{\mathrm{LAPs}}}=\\ & {\displaystyle \frac{\mathrm{1}}{m}}\sum _{i=\mathrm{1}}^m\left({\displaystyle \frac{(b\cdot {\mathrm{RF}}_{\mathrm{MODIS},\mathrm{anomaly}}^{\mathrm{LAPs}}(I_{\mathrm{LAPs}}{)}_i{)}^{\mathrm{2}}}{\begin{array}{c}(b\cdot {\mathrm{RF}}_{\mathrm{MODIS},\mathrm{anomaly}}^{\mathrm{LAPs}}(I_{\mathrm{LAPs}}{)}_i{)}^{\mathrm{2}}+(c\cdot {\mathrm{RF}}_{\mathrm{MODIS},\mathrm{anomaly}}^{\mathrm{LAPs}}(R_{\mathrm{eff}}{)}_i{)}^{\mathrm{2}}\\ +(d\cdot {\mathrm{RF}}_{\mathrm{MODIS},\mathrm{anomaly}}^{\mathrm{LAPs}}(\mathit{\theta }{)}_i{)}^{\mathrm{2}}\end{array}}}\right),\end{array}\end{array}$$

where m is the length of the data series. Similarly, we can obtain ${\mathrm{FC}}_{R_{\mathrm{eff}}}$ and FC_θ.

3.2.7 Calculation of in situ radiative forcing by LAPs in snow

${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ should be validated with measurements. However, due to the lack of radiative forcing measurements in NEC, we estimate the in situ radiative forcing (${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}})$ from measured BC_equiv values. Briefly, we use SNICAR to calculate the in situ reduction in snow albedo from BC_equiv and MODIS-retrieved R_eff. Then, the SBDART model is used to estimate ${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$.

4 Results

In January–February, seasonal snow widely covers northeastern China. For example, the area with a snow cover fraction of >50 % and snow duration period of >30 d is ∼75 % and ∼85 %, respectively (Fig. S5a and b), which is consistent with previous studies based on meteorological station data (Zhong et al., 2010) and satellite remote sensing data (Che et al., 2008). In addition, the area with a snow depth of >5 cm and snow water equivalent of >20 mm is ∼70 % and ∼70 %, respectively (Fig. S5c and d).

4.1 The spatial distribution of AOD and BC emissions

Northeastern China usually suffers from heavy local pollutant emissions, with high aerosol mass concentrations in winter (Wiedensohler et al., 2009). Figure 2a shows the spatial distribution of AOD at 550 nm derived from MODIS in NEC. We can find that AOD in the study areas ranges from 0.08 to 0.65 and shows strong spatial inhomogeneity. The largest AOD values are found in industrial areas in south-central NEC, where the largest urban areas in NEC are located with the major cities of Harbin, Changchun, and Shenyang. These areas are associated with the largest pollution emissions and anthropogenic activities in NEC (Wang et al., 2017). By comparison, the MODIS Aqua results show that the AOD in the west of NEC along the border of China is the smallest. Similar patterns of AOD were also found by Zhang et al. (2013) and Zhao et al. (2014). Previous studies indicated that BC represents the primary light-absorbing particles in the atmosphere (Cao et al., 2006) and seasonal snow (Wang et al., 2013). Figure 2b shows the spatial distribution of BC emission density during 2003–2014 in NEC. The pattern of BC emission density is very comparable to AOD, with the highest values of $>\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ g km⁻² month⁻¹ in south-central NEC and the lowest values of $<\mathrm{5}\times {\mathrm{10}}^{\mathrm{2}}$ g km⁻² month⁻¹ in the remote areas of northwestern China. Both the results of AOD and BC emission density imply that the seasonal snow in south-central NEC suffers from abundant BC deposition, and the radiative forcing by LAPs in snow is likely to be highest in NEC.

Figure 2 Spatial distribution of (a) MODIS AOD at 550 nm and (b) BC emission density in January–February in NEC. AOD data are from 2003 to 2017, and BC emission density data are from the research group at Peking University (http://inventory.pku.edu.cn/home.html, last access: 15 June 2018) from 2003 to 2014. The major cities in NEC are also shown in this figure.

4.2 The spatial distribution of snowfall frequency and land cover types

Snowfall is spatially varied in NEC and has a dominant effect on local fractional snow cover in the defined snow-covered areas, where we retrieved the radiative forcing by LAPs in snow in our study. Figure 3a shows the normalized snowfall frequency in January–February from 2003 to 2017. We can find that the highest snowfall frequency occurred in northwestern and southeastern NEC, where there are forest-covered areas (Fig. 3b). In contrast, the areas from central to southwestern NEC present the lowest snowfall frequency, which means that the fractional snow cover in these areas is likely to be lower than other areas and unable to reach the critical value that we used to define the snow-covered areas. On the other hand, land cover types will also affect the local fractional snow cover. From Fig. 3b, we can find that NEC presents spatially different land cover types; the mainland cover types are grasslands, croplands, and evergreen needleleaf (forests). Grasslands and croplands are mainly located in southwestern NEC and central NEC, respectively, while forests are distributed in northern and southeastern NEC. In our study periods, grasslands and croplands have limited influence on snow cover. However, in forest areas, even those completely covered by deep snow, forest will affect the derived surface reflectance from the MODIS Aqua satellite for the determination of snow-covered areas (further discussion in Sect. 5).

Figure 3 Spatial distribution of (a) the normalized snowfall frequency in January–February from 2003 to 2017 and (b) the different land cover types based on MODIS data in NEC. Snowfall data are from the ERA-Interim reanalysis. The major cities in NEC are also shown in this figure.

4.3 Radiative forcing by LAPs in snow

Figure 4 shows the identified snow-covered areas, which are primarily concentrated between 40 and 50^∘ N. Consistent with our analysis above, the low-snow-frequency areas of south-central and southwestern NEC and forest-covered areas of northern and southeastern NEC are not identified as snow-covered areas. According to the geographical distribution (Fig. 4a), we separated the studied area into three regions: western NEC (WNEC), central NEC (CNEC), and eastern NEC (ENEC).

Figure 4 The spatial distributions of average (a) I_LAPs, (b) R_eff, and (c) ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in NEC in January–February for 2003–2017. The background shows the spatial distribution of MODIS AOD values. The dotted areas are covered by forests. The major cities in NEC are also shown in this figure. According to the geographical distribution, we separate the study area into three regions: western NEC (WNEC), central NEC (CNEC), and eastern NEC (ENEC).

The spatial distributions of I_LAPs, R_eff, and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ are displayed in Fig. 4, and their statistics are presented in Fig. 5. On average, I_LAPs is $\sim \mathrm{0.27}\pm \mathrm{0.045}$, R_eff is $\sim \mathrm{261}\pm \mathrm{32}$ µm, and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is $\sim \mathrm{45.1}\pm \mathrm{6.8}$ W m⁻² in NEC. Regionally, ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is the largest and shows an average of $\sim \mathrm{50.9}\pm \mathrm{4.2}$ W m⁻² in CNEC, where industrial areas are located close to the largest urban areas in NEC; it therefore suffers from the most serious pollutant emissions among these three regions. ENEC displays the second largest radiative forcing with an average ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ of $\sim \mathrm{45.7}\pm \mathrm{4.5}$ W m⁻². The lowest value of $\sim \mathrm{41.0}\pm \mathrm{5.9}$ W m⁻² occurs in WNEC, where both AOD and BC emission densities are the lowest compared with the other two regions, which is not only due to the low local pollutant emissions, but also because the regional transport into this region in our study period is mostly from the clean northwest and suffers from little influence of human activities (Wang et al., 2015). For the individual regions, ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ presents an increase from north to south in CNEC that ranges from 40.4 to 64.6 W m⁻². In ENEC an east–west gradient of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is noted that ranges from 35.0 to 62.0 W m⁻². The most distinct intra-regional difference is in WNEC, where ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ ranges from 22.3 to 55.5 W m⁻². Generally, the patterns are consistent with those of AOD and BC emission density in NEC. Moreover, the spatial pattern of radiative forcing by LAPs in snow in this study is comparable with the results by Zhao et al. (2014), who first estimated the radiative forcing of LAPs in snow through the Weather Research and Forecasting (WRF) model and found that the radiative forcing in industrial source regions such as southern CNEC is obviously much higher than in border regions such as WNEC, which primarily resulted from the spatial differences of LAP dry and wet deposition. Compared with the results from other studies, Seidel et al. (2016) reported a radiative forcing of ∼20 W m⁻² in the Sierra Nevada in late February, which is lower than the result in NEC, even though the surface solar irradiance in the Sierra Nevada is higher. Painter et al. (2013b) reported an average radiative forcing of 215±63 W m⁻² in the Senator Beck Basin Study Area (SBBSA), SW Colorado, USA, which is approximately 4 times our retrieved radiative forcing near industrial areas in NEC. However, the snow grain size and the surface solar irradiance in their study period are larger than in our study by a factor of >2.5 and >4, respectively. These results indicate an abundant LAP content in the snow of CNEC. The regional and intra-regional patterns of variability in I_LAPs are quite similar to those of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$, which indicates the significant role of LAP content in determining the spatial distribution of radiative forcing; the average values of I_LAPs are $\sim \mathrm{0.311}\pm \mathrm{0.024}$ in CNEC, $\sim \mathrm{0.307}\pm \mathrm{0.026}$ in ENEC, and $\sim \mathrm{0.238}\pm \mathrm{0.031}$ in WNEC. In contrast to I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$, R_eff displays a smaller spatial variance and presents average values of $\sim \mathrm{285}\pm \mathrm{16}$, $\sim \mathrm{281}\pm \mathrm{15}$, and $\sim \mathrm{239}\pm \mathrm{29}$ µm in CNEC, ENCE, and WNEC, respectively. R_eff in WNEC is a little smaller compared with those in the other two regions, which is probably due to the higher snowfall frequency because a higher snowfall frequency indicates a longer duration of fresh finer snow at the surface (Wang et al., 2013; Seidel et al., 2016).

Figure 5 Statistics of average ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$, I_LAPs, and R_eff in NEC in January–February from 2003 to 2017.

Figure 6 shows the average uncertainties of I_LAPs, R_eff, and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ due to atmospheric correction in NEC in January–February from 2003 to 2017. The positive (negative) uncertainties of retrieved I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ from atmospheric correction are comparable and range from 9 % to 43 % (−10 % to −47 %) and 14 % to 57 % (−14 % to −47 %), respectively. Both I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ show larger uncertainties, as their values are smaller; the positive (negative) uncertainties of I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ are largest in WNEC and show averages of 21 % (−24 %) and 30 % (−28 %), while the lowest uncertainties of 13 % (−15 %) and 20 % (−20 %) for I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ are found in CNEC. It is because the uncertainty of snow albedo from atmospheric correction is almost similar in our study areas to that across the whole NEC region, as discussed in Sect. 3.6; however, the snow-albedo reduction is smaller in clean snow and larger in polluted snow, which results in a larger relative uncertainty of snow-albedo reduction in clean snow and a smaller relative uncertainty in polluted snow according to Eqs. (7) and (8). The positive (negative) uncertainties of R_eff are smaller compared with I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$, and they range from 14 % to 18 % (−12 % to −16 %), which is comparable with the errors between MODIS-retrieved and in situ measured snow grain size discussed in Sect. 3.2.2. Moreover, the uncertainties are spatially quite consistent for R_eff, which is different from I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. We note that the positive and negative uncertainties of all I_LAPs, R_eff, and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ are asymmetric, which is primarily due to the nonlinear characteristics of radiative transfer in the SNICAR model (Painter et al., 2007).

Figure 6 (a) Negative and (b) positive uncertainty of average I_LAPs in NEC in January–February from 2003 to 2017. Panels (c) and (d) are similar to (a) and (b), but for R_eff. Panels (e) and (f) are similar to (a) and (b), but for ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. The background shows the spatial distribution of MODIS AOD values. The dotted areas are covered by forests. The major cities in NEC are also shown in this figure.

As discussed in Sect. 3, the spatial distribution of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ depends on LAPs, R_eff, and θ. Previous studies have attempted to retrieve the radiative forcing by LAPs in snow by using remote sensing (e.g., Painter et al., 2012a, 2013b); however, attributions of the spatial variations of radiative forcing by LAPs in snow are sparse and need to be further investigated. Therefore, we would like to qualify the contribution of each factor to the spatial variance of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. Combing sensitivity tests and the method of Huang and Yi (1991) as discussed in Sect. 3.2.6, we estimate the fractional contribution of I_LAPs (the indicator of LAPs), R_eff, and θ to the spatial variance of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in our study areas across NEC (Fig. 7). We can find that the contribution from LAPs is largest with a value of 74.6 %, while R_eff and θ make contributions of 21.2 % and 4.2 %, respectively, in NEC. The result indicates that the LAP content in snow plays a dominant role in determining the spatial distribution of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. Regionally, the contribution of LAPs in WNEC (62.1 %) is smaller than those of 73.9 % and 83.4 % in CNEC and ENEC, while R_eff shows a higher contribution of 28.1 % in WNEC than those of 19.6 % and 13.9 % in CNEC and ENEC. The results point out a less important effect of LAPs but a more important effect of R_eff on the spatial distribution of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in WNEC compared with those in CNEC and ENEC. In addition, the contribution of θ is smaller in ENCE (2.7 %) than those of 9.8 % and 6.5 % in WNEC and CNEC, which is primarily due to the smallest altitude range of ENEC among those three regions.

Figure 7 Fractional contribution of average I_LAPs (the indicator of LAPs), R_eff, and solar zenith angle (θ) to the spatial variance of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in January–February for 2003–2017.

Seidel et al. (2016) reported that the variations in LAP contents in snow are dominated by LAP deposition and snowfall. Previous studies have also reported that BC is the dominant LAP type in NEC (Wang et al., 2013). Zhao et al. (2014) simulated LAP contents and their radiative forcing in seasonal snow using WRF-Chem coupled with the SNICAR model and indicated that the radiative forcing by LAPs in snow in NEC is primarily due to BC. Therefore, to examine the spatial distributions of retrieved I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$, we display the distribution of snowfall (Fig. 3a) and BC dry and wet deposition (Fig. 8). BC dry deposition is highest in the largest urban areas of NEC with the major cities of Harbin, Changchun, and Shenyang, then decreases sharply outwards from the central urban areas to remote areas (Fig. 8a). Different from BC dry deposition, which is dominated by BC concentrations in the atmosphere, BC wet deposition is affected by both BC concentrations and precipitation and shows an increase from northwest to southeast. Generally, the spatial patterns of BC dry and wet deposition are similar to I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. For example, areas with higher BC dry and wet deposition, such as industrial polluted NEC, show higher I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. Moreover, from Fig. 9a–c, we can find that the correlations of I_LAPs with BC dry and wet deposition and snowfall (R²=0.81, 0.73, and 0.14) are all significant at the 99 % confidence level. The correlations of I_LAPs with BC dry and wet deposition in WNEC are relatively lower than those in CNCE and ENEC, which is partly due to the effect of dust in this region (Wang et al., 2013; Zhao et al., 2014). Furthermore, using the method of multiple linear regression, we fitted I_LAPs using BC dry and wet deposition and snowfall data. Figure 9d shows the scatterplots of I_LAPs and fitted I_{LAPs_fit}. We can find that I_{LAPs_fit} is highly correlated with I_LAPs, and BC dry and wet deposition and snowfall could explain 84 % of the spatial variance of I_LAPs. The result confirms the reasonability of the spatial patterns of retrieved I_LAPs and thus ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in NEC. In addition to MERRA-2 BC deposition data and ERA-Interim snowfall data used in Fig. 9, we also used other types of BC deposition and snowfall data to fit I_LAPs. Table S1 in the Supplement shows the R² between MODIS-retrieved I_LAPs and fitted I_{LAPs_fit} based on different datasets as discussed in Sect. 2.3 and 2.4. The values of R² are very similar in a range of 0.81–0.84, which further indicates that the spatial pattern of retrieved I_LAPs is reasonable and independent of the data types used for validation.

Figure 8 Spatial distribution of average (a) dry and (b) wet deposition of BC in NEC in January–February from 2003 to 2017. BC deposition data are from the MERRA-2 reanalysis.

Figure 9 Scatterplots of I_LAPs versus (a) BC dry deposition, (b) BC wet deposition, (c) normalized snowfall frequency, and (d) fitted I_LAPs (I_{LAPs_fit}), which is fitted with BC dry and wet deposition and snowfall frequency using multiple linear regression. BC deposition data are from the MERRA-2 reanalysis, and snowfall data are from the ERA-Interim reanalysis in January–February from 2003 to 2017.

4.4 Comparisons of MODIS-retrieved and in situ estimated radiative forcing by LAPs in snow

Figure 10 shows the distribution of the sample sites and the measured BC_equiv concentration in surface snow at each site. Circles and squares represent the snow samples collected in 2010 (Wang et al., 2013) and 2014 (Wang et al., 2017), respectively. Generally, BC_equiv concentration ranges mostly from ∼0.1 to ∼3.0 µg g⁻¹ and shows an increase from northwest to southeast. The highest BC_equiv concentration is found in CNEC, while the lowest is in WNEC. Figure 11a displays a comparison of MODIS-retrieved radiative forcing (${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$) and in situ radiative forcing (${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}})$ estimated based on measured BC_equiv concentration. In general, the mean absolute error (MAE) for ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ against ${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$ is 15.3 W m⁻². The ratios of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ to ${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$ ($R_{\mathrm{in}\mathrm{situ}}^{\mathrm{MODIS}}$) fall mainly in the range of 1–2. The errors indicate larger positive biases at lower ${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$ values, whereas smaller biases are noted at higher ${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$ values. The results of this bias analysis are comparable with those reported by Painter et al. (2012a). Figure 11b shows a scatterplot of $R_{\mathrm{in}\mathrm{situ}}^{\mathrm{MODIS}}$ versus BC_equiv. We can find that $R_{\mathrm{in}\mathrm{situ}}^{\mathrm{MODIS}}$ and BC_equiv display a good correlation; the best-fitting equation is $R_{\mathrm{in}\mathrm{situ}}^{\mathrm{MODIS}}=\mathrm{1.690}\cdot {\mathrm{BC}}_{\mathrm{equiv}}^{-\mathrm{0.522}}$, and the R² is 0.89 (99 % confidence level). This result indicates that the biases in the ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ retrievals are negatively correlated with the LAP concentrations in NEC. Considering that the typical concentration of BC_equiv in clean snow in NEC is 0.15 µg g⁻¹, the bias in ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ can be as high as 350 % in some areas, such as WNEC. In other areas with very polluted snow, such as southern CNEC (where the BC_equiv values are typically 2.5 µg g⁻¹), the bias is ∼5 %. Thus, considering the values reported by Wang et al. (2013, 2017), the biases in ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ largely fall in the range of ∼5 % to ∼350 % in NEC. Comparing Fig. 11 with Fig. 6, we find that the biases in ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in polluted snow are comparable with the uncertainties of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ due to atmospheric corrections. However, in clean snow, the uncertainties from atmospheric corrections could not sufficiently explain the biases in retrieved ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. There are three probable reasons. (a) For clean snow, retrieved radiative forcing is very sensitive to MODIS-derived surface snow reflectance (Eq. 4); although we have corrected the errors of snow reflectance from the protrusion of vegetation in our study areas of high snow cover fractions, the uncertainties from fractional snow cover (FSC) calculation and the vegetation reflectance will effectively influence the corrections of snow reflectance (Eq. 5). (b) Painter et al. (2012b) validated the retrieved radiative forcing by LAPs in snow in the upper Colorado River Basin using in situ estimates based on radiation towers, and they also found that the biases in the case of low radiative forcing could be up to severalfold. They pointed out that MODIS cannot proceed with a continuous spectral measurement of a continuous variable forcing like what LAPs afford to snow albedo due to the variably spaced and discrete bands of MODIS, which prevents a more quantitative retrieval and thus results in a non-negligible uncertainty in radiative forcing retrieval. (c) We use the average of MODIS-retrieved radiative forcing in a pixel size of 0.05^∘ × 0.05^∘ to compare with the in situ radiative forcing calculated using the observed BC_equiv concentration with the sample site located in the center of the pixel. Such a comparison may not be true at some sites due to the inhomogeneous spatial distribution of snow and LAP contents, which will influence radiative forcing estimates, especially in clean snow (Zhao et al., 2014). Therefore, we note that the number of sample sites is still limited and more field campaigns are needed to validate the accuracy of MODIS retrievals and then correct the retrieved radiative forcing.

Figure 10 Spatial distribution of the measured BC_equiv concentration in surface snow in NEC. Circles and squares represent the snow samples collected in 2010 (Wang et al., 2013) and 2014 (Wang et al., 2017), respectively.

4.5 Limitations

The determination of snow-covered areas represents a limitation of the method used in this study, which restricts our study to areas with high snow cover fractions; thus, we cannot estimate ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ across NEC as a whole. In the future, we will attempt to apply other satellite data with higher spatial resolution and use the spectral differences between different land cover types to distinguish the spectral reflectance of snow in mixed pixels. These improvements will permit us to expand our work to areas with limited snow cover. Another limitation is that we retrieve only the instantaneous radiative forcing at the surface under clear-sky conditions at the time of MODIS overpass, and these measurements do not represent a time-integrated average over the studied period. However, the estimation of temporally resolved radiative forcing is much more difficult given the significant effects of clouds, atmospheric components, θ, and the time-varying snow reflectance.

Figure 11 Scatterplots of (a) ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ versus ${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$ and (b) $R_{\mathrm{in}\mathrm{situ}}^{\mathrm{MODIS}}$ versus BC_equiv.

5 Discussion

In our study, we did not retrieve the radiative forcing in the northern and southeastern parts of NEC. In those regions, snowfall is frequent, the percent of snow cover is very high, and snow is also very deep. For example, in northern NEC, the averaged snow depth is ∼20 cm, and in the areas near the Changbai Mountains of southeastern NEC, snow depth could be up to ∼40 cm (Wang et al., 2013). However, due to the presence of forest cover, the reflected radiation received by sensors aboard satellites in those areas is mostly due to trees. For example, Fig. 12 shows the true color map of MODIS in NEC on 23 January 2010, and we can see that in the northern and southeastern parts of NEC, the observed objects from MODIS are almost all trees, not the snowpack under trees, although snow is almost completely covering the area (Wang et al., 2013). Therefore, in those forest areas, discussing the radiative forcing by LAPs in snow is extremely difficult due to the influence of trees. Bond et al. (2006) also indicated that LAPs in snow masked by forests contribute little to radiative forcing. They further pointed out that model representation of and forcing sensitivity to cover ranges of forests have not been verified, and this is a boundless uncertainty in modeling radiative forcing by LAPs in snow at present. However, most modeling studies that have simulated the radiative forcing by LAPs in snow did not take trees into consideration and estimated the radiative forcing over entire boreal forest areas in the Northern Hemisphere. For example, Flanner et al. (2007) applied the SNICAR model coupled a general circulation model to estimate the radiative forcing and response from BC in snow-covered areas over the whole Northern Hemisphere. Nevertheless, due to the presence of trees in the extensive boreal forest areas, the simulated radiative forcing is unreal as the incident radiation is reflected by trees but not by the snowpack. Zhao et al. (2014) simulated BC and dust and their radiative forcing in seasonal snow in northern China. They found that the radiative forcing by BC and dust is very high in southeastern NEC, where there are forest areas. But, in fact, in those areas the simulated radiative forcing by LAPs is also unreal. Therefore, we note that estimating the radiative forcing by LAPs in forest areas should consider the influence of trees.

Figure 12 A true color map of MODIS in NEC on 23 January 2010.

6 Conclusions

In this study, we retrieve I_LAPs, R_eff, and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ across NEC in January–February from 2003 to 2017 using MODIS data, together with a snow-albedo model (SNICAR) and a radiative transfer model (SBDART). On average, I_LAPs is $\sim \mathrm{0.27}\pm \mathrm{0.045}$, R_eff is $\sim \mathrm{261}\pm \mathrm{32}$ µm, and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ is $\sim \mathrm{45.1}\pm \mathrm{6.8}$ W m⁻² in NEC. The distribution of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ presents distinct spatial differences; the lowest value is 22.3 W m⁻², which occurs in remote western NEC, and the highest value is 64.6 W m⁻², which occurs near the industrial areas in central NEC. Both I_LAPs and ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ show larger uncertainties from atmospheric correction, as their values are smaller. We make a first attempt to attribute the variations of radiative forcing based on remote sensing. The results point out that I_LAPs, R_eff, and θ make fractional contributions of 74.6 %, 21.2 %, and 4.2 % to the spatial variance of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in our study areas across NEC. The result confirms that the LAP content in snow plays a dominant role in determining the spatial distribution of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$. We also analyze the distribution of BC dry and wet deposition and snowfall, and we find that they could explain 84 % of the spatial variance of I_LAPs, which indicates the reasonability of the spatial patterns of I_LAPs and thus ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ in NEC. Finally, we validate the retrieved ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ values using in situ estimated radiative forcing (${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$). The mean absolute error (MAE) of ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ against ${\mathrm{RF}}_{\mathrm{in}\mathrm{situ}}^{\mathrm{estimated}}$ is 15.3 W m⁻². The biases in the ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ retrievals display a negative correlation with the LAP concentrations in NEC. Considering typical concentrations of BC_equiv, which range from ∼0.15 to ∼2.5 µg g⁻¹, the biases in ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$ fall primarily within the range of ∼5 % to ∼350 % in NEC.