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.1 Models

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.

Download

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 2Spatial 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 3Spatial 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 4The 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 5Statistics of average ${\mathrm{RF}}_{\mathrm{MODIS}}^{\mathrm{LAPs}}$, I_LAPs, and R_eff in NEC in January–February from 2003 to 2017.

Download

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

Download

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 8Spatial 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 9Scatterplots 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.

Download

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 10Spatial 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 11Scatterplots 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.

Download