Simulating black carbon and dust and their radiative forcing in seasonal snow : a case study over North China with field campaign measurements

A state-of-the-art regional model, the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008) coupled with a chemistry component (Chem) (Grell et al., 2005), is coupled with the snow, ice, and aerosol radiative (SNICAR) model that includes the most sophisticated representation of snow metamorphism processes available for climate study. The coupled model is used to simulate black carbon (BC) and dust concentrations and their radiative forcing in seasonal snow over North China in January–February of 2010, with extensive field measurements used to evaluate the model performance. In general, the model simulated spatial variability of BC and dust mass concentrations in the top snow layer (hereafter BCS and DSTS, respectively) are consistent with observations. The model generally moderately underestimates BCS in the clean regions but significantly overestimates BCS in some polluted regions. Most model results fall within the uncertainty ranges of observations. The simulated BCS and DSTS are highest with > 5000 ng g −1 and up to 5 mg g−1, respectively, over the source regions and reduce to < 50 ng g −1 and < 1 μg g−1, respectively, in the remote regions. BCS and DSTS introduce a similar magnitude of radiative warming ( ∼ 10 W m−2) in the snowpack, which is comparable to the magnitude of surface radiative cooling due to BC and dust in the atmosphere. This study represents an effort in using a regional modeling framework to simulate BC and dust and their direct radiative forcing in snowpack. Although a variety of observational data sets have been used to attribute model biases, some uncertainties in the results remain, which highlights the need for more observations, particularly concurrent measurements of atmospheric and snow aerosols and the deposition fluxes of aerosols, in future campaigns.


Introduction
Snow is an important component in the Earth's climate system.Water from snowmelt generates runoff to fill rivers and reservoirs in many regions of the world.Snow cover can change the surface albedo of grassland by a factor of 3-4 and forested regions by a factor of 2-3 (Thomas and Rowntree, 1992;Betts and Ball, 1997).Therefore, snow, at a larger scale, regulates the temperature of the Earth's surface and alters the general circulation, and, at a smaller scale, affects regional climate (e.g., Barnett et al., 1988;Walland and Simmonds, 1996).Small changes in snow albedo can have large impacts on surface warming due to the rapid feedbacks involving changes to snow morphology, sublimation, and melt rates (Bond et al., 2013).Snow albedo can be influenced by many factors including snow depth, snow grain size, and snow impurity.Snow impurity is mainly due to the light absorbing aerosols (LAA) deposited on/in the snowpack, Published by Copernicus Publications on behalf of the European Geosciences Union.
including primarily black carbon (BC), brown carbon (BrC), and dust (Warren andWiscombe, 1980, 1985;Hansen and Nazarenko, 2004;Painter et al., 2007;Jacobson, 2004Jacobson, , 2010;;Flanner et al., 2007Flanner et al., , 2009Flanner et al., , 2012;;Skeie et al., 2011;Lin et al., 2014;Wang et al., 2014).A recent estimate for global effective BC forcing in snow and ice is +0.13 W m −2 (Bond et al., 2013).Dust residing in snow can also reduce snow albedo, particularly for the visible bands.In the upper Colorado River basin, a total end-of-year dust concentration of 0.23-4.16mg g −1 in the top 3 cm of the snow column can result in mean springtime dust radiative forcing in snow from 31 to 49 at the alpine site and 45 to 75 W m −2 at the subalpine site (Painter et al., 2012;Skiles et al., 2012).
The LAA effects in snow could significantly affect the atmospheric general circulation and hydrological cycle at global and regional scales in multiple ways (Qian et al. 2009(Qian et al. , 2011;;Painter et al., 2010;Skiles et al., 2012;Bond et al, 2013;Sand et al., 2013).For example, over Asia, BC in snow could increase the surface air temperature and reduce spring snowpack over the Tibetan Plateau (TP) (Qian et al., 2011).Earlier snowmelt can shift the timing of peak runoff and its duration and modulate the soil moisture and surface fluxes.These impacts may result in an earlier onset of the South Asian Monsoon and an increase of moisture, cloudiness, and convective precipitation over northern India (Qian et al., 2011).They may also affect the East Asian monsoon, resulting in changes in summer precipitation pattern (Qian et al., 2011).In the western USA, Qian et al. (2009) found that BC deposition on the mountain snowpack could also cause a reduction in snow accumulation and an earlier snowmelt that impacts water resources during the drier summer months.In the upper Colorado River basin, dust radiative forcing in snow could shorten the snow cover duration by 21 to 51 days, cause the peak runoff to occur 3 weeks earlier on average, and reduce the annual runoff by more than ∼ 5 % (Painter et al., 2010;Skiles et al., 2012).
Although there have been some modeling studies estimating the LAA concentrations in snow and their impact on snow albedo, solar absorption in snow, and climate and hydrological cycle, most previous studies used global climate models (e.g., Jacobson, 2004Jacobson, , 2010;;Flanner et al., 2007Flanner et al., , 2012;;Skeie et al., 2011;Qian et al., 2011;Bauer and Menon, 2012;Lin et al., 2014) or snow energy balance models coupled with simplified radiation schemes with observed LAA concentrations in snow as inputs (e.g., Painter et al., 2010Painter et al., , 2012;;Skiles et al., 2012).Both amounts and properties of snow and aerosols are extremely heterogeneous, particularly in regions with complex topography (e.g., Painter et al., 2010;Zhao et al., 2013a).Global climate models with relatively coarse horizontal resolutions (typically 1-2 • ) cannot fully resolve the spatial and temporal variability of aerosols in snow.In addition, global climate models with some simplified treatments of physics and chemistry processes may introduce additional uncertainties in simulating snow and aerosols (Zhao et al., 2013b, c).
Few studies of the climatic impacts of LAA in snowpack used regional climate models.Although regional models are more computationally demanding, the ability to resolve complex terrain using high spatial resolution is particularly important for simulating mountain precipitation and snowpack (e.g., Leung et al., 2003;Leung and Qian, 2003;Rasmussen et al., 2011;Yoon et al., 2012).More sophisticated physics and chemistry treatments used in some regional models can provide additional benefits.Although Qian et al. (2009) used a regional modeling framework to investigate the impact of BC in snow over the western USA, they estimated the BC content in snow based on a single year simulation with a coupled chemistry component and applied the same BC deposition to simulations of multiple years to estimate BC-insnow effects.This method limits the interactions between BC deposition and climate that could influence the BC-induced snow albedo perturbations.Furthermore, snow albedo perturbations from snow aging were ignored, which limits the potential for accelerated snow aging induced by a BC-in-snow effect.Lastly, they did not simulate the dust-in-snow effect, which likely had a much greater impact than BC in some regions (e.g., remote mountains) in specific seasons (Painter et al., 2007).
Although SNICAR was released with the CLMv4.0 in WRF-Chem (v3.5.1), it is not connected with the atmospheric aerosol deposition and, therefore, the aerosol radiative effect in snowpack does not function.This study couples the WRF-Chem simulated aerosols with SNICAR and makes an aerosol radiative effect in snowpack function.The modification of SNICAR by Flanner et al. (2012) is also implemented into the model.The coupled model is used to simulate the BC and dust concentrations and their radiative forcing in snowpack during a previous field campaign that collected seasonal snow on a road trip at 46 sites over North China in January-February of 2010 (Huang et al., 2011).Brown carbon is excluded from the SNICAR simulations in this study because its optical properties are poorly constrained (Flanner et al., 2007(Flanner et al., , 2009)).This study represents the first effort on evaluating WRF-Chem for simulating BC and dust in snowpack at a relatively high spatial resolution against the field campaign measurements of BC content in snow, and on estimating the radiative forcing of BC and dust in snow over North China.The rest of the paper is organized as follows.Sections 2 and 3 provide detailed description of the WRF-Chem and SNICAR models and the observations used in this study.In Sect.4, the simulated spatial and temporal variation of snow and aerosol content in snow are evaluated against measurements and characterized.Then, the radiative forcing of BC and dust in snow over North China are estimated.The findings are concluded in Sect. 5.

Model
In this study, WRF-Chem (v3.5.1), briefly described in Sect.2.1, is used to couple the SNICAR model with interactive aerosols as described in Sect.2.2 below.Section 2.3 discusses the setup of model simulations.The emissions used in the simulations including anthropogenic and biomass burning emissions and online calculated emissions of mineral dust and sea salt are described in Sect.2.4.

WRF-Chem
In this study, the MOSAIC (Model for Simulating Aerosol Interactions and Chemistry) aerosol model (Zaveri et al., 2008) coupled with the CBM-Z (carbon bond mechanism) photochemical mechanism (Zaveri and Peters, 1999), one of the chemistry options in WRF-Chem, is used and coupled with the SNICAR model.The MOSAIC aerosol scheme uses the sectional approach to represent aerosol size distributions with a number of discrete size bins, either four or eight bins in the current version of WRF-Chem (Fast et al., 2006).All major aerosol components including sulfate (SO 2− 4 ), nitrate (NO − 3 ), ammonium (NH + 4 ), black carbon (BC), organic matter (OM), sea salt, and mineral dust are simulated in the model.The MOSAIC aerosol scheme includes physical and chemical processes of nucleation, condensation, coagulation, aqueous phase chemistry, and water uptake by aerosols.Dry deposition of aerosol mass and number is simulated following the approach of Binkowski and Shankar (1995), which includes both particle diffusion and gravitational effects.Wet removal of aerosols by gridresolved stratiform clouds/precipitation includes in-cloud removal (rainout) and below-cloud removal (washout) by impaction and interception, following Easter et al. (2004) and Chapman et al. (2009).In this study, cloud-ice-borne aerosols are not explicitly treated in the model but the removal of aerosols by the droplet freezing process is considered.Convective transport and wet removal of aerosols by cumulus clouds follow Zhao et al. (2013c).
Aerosol optical properties such as extinction, singlescattering albedo (SSA), and asymmetry factor for scattering are computed as a function of wavelength for each model grid box.Aerosols are assumed internally mixed in each bin, i.e., a complex refractive index is calculated by volume averaging for each bin for each chemical constituent of aerosols.The Optical Properties of Aerosols and Clouds (OPAC) data set (Hess et al., 1998) is used for the shortwave (SW) and longwave (LW) refractive indices of aerosols, except that a constant value of 1.53+0.003i is used for the SW refractive index of dust following Zhao et al. (2010Zhao et al. ( , 2011)).A detailed description of the computation of aerosol optical properties in WRF-Chem can be found in Fast et al. (2006) and Barnard et al. (2010).Aerosol radiative feedback is coupled with the Rapid Radiative Transfer Model (RRTMG) (Mlawer et al., 1997;Iacono et al., 2000) for both SW and LW radiation as implemented by Zhao et al. (2011).The optical properties and direct radiative forcing of individual aerosol species in the atmosphere are diagnosed following the methodology described in Zhao et al. (2013a).Aerosol-cloud interactions were included in the model by Gustafson et al. (2007) for calculating the activation and re-suspension between dry aerosols and cloud droplets.

SNICAR
The SNICAR is a multilayer model that accounts for vertically heterogeneous snow properties and heating and influence of the ground underlying snow (Flanner and Zender, 2005;Flanner et al., 2007Flanner et al., , 2009Flanner et al., , and 2012)).The model uses the theory from Wiscombe and Warren (1980) and the two-stream, multilayer radiative approximation of Toon et al. (1989).SNICAR can well simulate snow surface albedo as well as the radiative absorption within each snow layer.It can also simulate aerosol content and radiative effect in snow, and was first used to study the aerosol heating and snow aging in a global climate model by Flanner et al. (2007).The SNICAR simulated change of snow albedo for a given BC concentration in snow has been validated with recent laboratory and field measurements (Brandt et al., 2011;Hadley and Kirchstetter, 2012).More detailed description of the SNICAR model can be found in Flanner and Zender (2005) and Flanner et al. (2007Flanner et al. ( , 2012)).In WRF-Chem with CLM chosen as the land surface scheme, SNICAR defines radiative layers that match the five thermal layers in CLM that vertically resolve snow thermal processes, densification, and meltwater transport (e.g., Oleson et al., 2010).Here, we focus on the description of processes associated with modeling aerosols and their optical properties in snowpack.

Aerosol deposition and scavenging
BC and dust can deposit on/in snow through dry and wet deposition processes as discussed above.In this study, atmospheric aerosol wet deposition is calculated with the prognostic precipitation, thus avoiding inconsistencies associated with prescribed aerosol deposition and meteorology.Following Flanner et al. (2012), BC can mix externally and internally with precipitation hydrometeors while dust can only mix externally with precipitation hydrometeors.BC and dust removed through dry deposition processes (sedimentation and turbulent mix-out), often the dominant sink of large particles near emission sources, are mixed externally with snow grains.Dust removed through wet deposition is also treated as being externally mixed with snow grains.BC removed through below-cloud scavenging, where particles are collected by falling hydrometeors through Brownian motion, electrostatic forces, collision, and/or impaction, is also assumed to remain attached to the outside of the hydrometeors rather than becoming internalized.Cloud-borne BC collected by stratiform precipitation is treated as being internally mixed with snow grains.Interstitial and cloud-borne BC removed through convective scavenging processes is assumed to mix externally and internally, respectively, with snow grains.Although determining the mixing state of BC deposited through convective scavenging is uncertain, only a few percent of precipitation falling on snow surfaces comes from convection (Flanner et al., 2007;Sect. 4.2.1).Following Flanner et al. (2012), we do not simulate the distribution of aerosol number concentrations within the precipitating hydrometeors, which may influence total absorption.
In this study, SNICAR in WRF-Chem simulates four tracers of dust and two tracers of BC in snow following Flanner et al. (2012).The four tracers of dust represent four sizes of dust with diameters of 0.1-1, 1-2.5, 2.5-5, and 5-10 µm.All four dust tracers are assumed to mix externally with snow grains.With the MOSAIC aerosol model, WRF-Chem can simulate dust in the atmosphere with four size bins (i.e., 0.039-0.156,0.156-0.625,0.625-2.5, and 2.5-10.0µm in dry diameter) or eight size bins (0.039-0.078, 0.078-0.156,0.156-0.312,0.312-0.625,0.625-1.25,1.25-2.5,2.5-5.0, and 5.0-10.0µm in dry diameter).Both of them are coupled with the dust tracers in SNICAR.The coupling of dust size distributions between WRF-Chem and SNICAR are summarized in Table 1.Although both four-and eight-bin size representations in MOSAIC are coupled with SNICAR, the analysis presented in this paper focuses on the simulations using the eight-bin size representation because it is more accurate in representing aerosol size distributions (Zhao et al., 2013c).The two BC tracers in SNICAR represent interstitial and within-ice BC, respectively, in order to account for the light absorption enhancement by snowpack containing BC particles residing within snow grains (will be discussed in Sect.2.2.2).Following Flanner et al. (2012), this study fixes two BC tracers in snow with 0.1 µm dry effective radius, al-though size distributions of BC may influence the BC absorption in snow, which deserves further investigation.
The scavenging of aerosols in snow by meltwater is simulated by SNICAR.CLM accounts for meltwater drainage by adding excess water to the layer beneath when the liquid content exceeds the layer's holding capacity, defined by snow porosity and irreducible water saturation.SNICAR assumes aerosol inclusion if melt water is proportional to its mass mixing ratio multiplied by a scavenging factor.Thus the change rate of aerosol mass in each layer i is where m i is the absolute mass of aerosol in layer i, k is the scavenging ratio, q i is the mass flux of water out of layer i, c i is the mass mixing ratio of aerosol in layer i (aerosol mass divided by liquid and solid water mass), and D is the sum of wet and dry deposition of atmospheric aerosols, added only to the surface layer of snowpack.After deposition, aerosols are mixed instantly and uniformly in the model surface layer, which never exceeds 3 cm thick.In this study, following Flanner et al. (2012), k for interstitial BC (k pho ) and within-ice BC (k phi ) are assumed to be 0.03 and 0.20, respectively, and k for four dust tracers are 0.02, 0.02, 0.01, and 0.01.These scavenging ratios are highly uncertain and need to be constrained by observations in the future (Flanner et al., 2012;Qian et al., 2014), but the values used here are generally reasonable compared to observations in the high latitudes (Doherty et al., 2013).In addition, as fresh snowfall is added, CLM continuously divides the model surface snow layer and prevents excessive aerosol accumulation.Therefore, the concentrations of each tracer within each snow layer depend on atmospheric deposition rates, new snowfall, ice sublimation, meltwater flushing, and layer combinations and divisions (Flanner et al., 2007(Flanner et al., , 2012;;Oleson et al., 2010).

Optical properties
Optical properties of BC mixed internally and externally with snow grains are treated separately in SNICAR (Flanner et al., 2012).In this study, SNICAR accounts for the light absorption by snowpack containing BC particles residing within snow grains, which could significantly increase the solar absorption compared to interstitial BC.The refractive indices of BC from Chang and Charalampopoulos (1990) and the effective radius of 0.1 µm are used for offline calculated Mie parameters for BC (i.e., extinction, SSA, and scattering asymmetry parameters) (Flanner et al., 2007 and2012).Building on Chylek et al. (1984), the dynamic effective medium approximation method is applied to derive absorption enhancement ratio over broad size ranges of snow grains and BC expected for snowpack.Then, the mass absorption cross section of BC internally mixed with snow grains is scaled by the derived enhancement ratio, which is determined as a function of BC and snow effective radius online via a look-up  table.Application of the enhancement ratio was justified as a reasonable approach to account for absorption enhancement by BC internally mixed with snow grains (Flanner et al., 2012).Dust optical properties in snowpack are obtained using the Maxwell-Garnett mixing approximation along with the Mie theory.Dust in snowpack is assumed to be a blend of SiO 2 (47.6 % by volume), limestone (2 %), montmorillonite (25 %), illite (25 %), and hematite (0.4 %) with four size bins following the lognormal distribution with a number median radius of 0.414 µm and standard deviation of 2 (Flanner et al., 2009).Since dust is radiatively active, its co-existence could decrease the absorption by BC in snowpack (Flanner et al., 2009).Accounting for enhanced absorption by dust mixed internally with snow grains could further decrease BC absorption, but this is not considered in this study following Flanner et al. (2012).SNICAR treats the snow particle as a collection of ice spheres, with effective number median radius of 30-1500 µm.Mie parameters for one visible (0.3-0.7 µm) and four near-infrared (NIR) (0.7-1.0, 1.0-1.2,1.2-1.5, and 1.5-5.0µm) wave bands are computed offline for lognormal distributions of this wide size range of snow grains for computational efficiency.For radiative transfer calculations, the snow grains and six aerosol tracers are treated as external mixtures (with absorption enhancements applied to withinice BC), by summing the extinction optical depths of each component, weighting the individual SSA by optical depths, and weighting the asymmetry parameters by the product of optical depths and SSA (Flanner et al., 2007(Flanner et al., , 2012)).A detailed description of the computation of optical properties of snow with aerosols in SNICAR can be found in Flanner et al. (2012).

Numerical experiments
In this study, the WRF-Chem simulation is performed at 36 km × 36 km horizontal resolution with 120 × 90 grid cells (95-140 • E, 30-60 • N) (Fig. 1) and 35 vertical layers up to 100 hPa.Although only the results for January-February of 2010 were analyzed, the simulation was conducted from 1 October 2009 to 25 February 2010 to allow accumulation of snowpack and aerosols in snowpack.The meteorological initial and lateral boundary conditions are derived from the National Center for Environmental Prediction final analysis (NCEP/FNL) data at 1 • horizontal resolution and 6 h temporal intervals.The modeled u component and v component wind and atmospheric temperature are nudged towards the NCEP/FNL reanalysis data with a nudging timescale of 6 h (Stauffer and Seaman, 1990).The MYJ (Mellor-Yamada-Janjic) planetary boundary layer scheme, CLM land surface scheme, Morrison 2-moment microphysics scheme, Kain-Fritsch cumulus scheme, and RRTMG longwave and shortwave radiation schemes are used in this study.The chemical initial and boundary conditions are provided by a quasiglobal WRF-Chem simulation for the same time period to include long-range transported chemical species.The quasiglobal WRF-Chem simulation is performed at 1 • × 1 • horizontal resolution using a quasi-global channel configuration (using periodic boundary conditions in the zonal direction) with 360 × 130 grid cells (180 The quasi-global simulation is configured in a way similar to the regional simulation and also driven by the NCEP/FNL data.More details about the quasi-global WRF-Chem simulation can be found in Zhao et al. (2013c).

Emissions
Anthropogenic emissions are obtained from the Asian emission inventory described by Zhang et al. (2009) (Ginoux et al., 2001), and the emitted dust particles are distributed into the MOSAIC aerosol size bins following a theoretical expression based on the physics of scale-invariant fragmentation of brittle materials derived by Kok (2011).For MOSAIC 8-bin, dust par-ticles are emitted into eight size bins with mass fractions of 10 −6 , 10 −4 , 0.02, 0.2, 1.5, 6, 26, and 45 %.For MO-SAIC 4-bin, dust particles are emitted into four size bins with mass fractions of 10 −4 , 0.22, 7.5, and 71 %, respectively.WRF-Chem has been widely used to simulate the dust life cycle and climatic impact at a global scale (Zhao et al., 2013a) and regional scales over West Africa (Zhao et al., 2010(Zhao et al., , 2011)), Saudi Arabia (Kalenderski et al., 2013), North America (Zhao et al., 2012), and East Asia (Chen et al., 2013(Chen et al., , 2014)).The BC and dust emissions within the model domain are indicated in Fig. 1.
For the quasi-global WRF-Chem simulation that provides the chemical boundary for the regional simulation, anthropogenic emissions for the year 2000 are obtained from the Reanalysis of the TROpospheric (RETRO) chemical composition inventories (ftp://ftp.retro.enes.org/pub/emissions/aggregated/anthro/0.5x0.5/)except over East Asia and the United States, where anthropogenic emissions are from the Asian 2006 emission inventory (Zhang et al., 2009) and from the US National Emission Inventory (NEI) 2005 (WRF-Chem user guide from http://ruc.noaa.gov/wrf/WG11/Users_guide.pdf),respectively.Emissions of biomass burning aerosols, sea salt, and dust are treated in the same way as described above for the regional simulation.

Snow and its BC content data set
The primary observational data set used in this study to evaluate the model simulations of snow and its BC content is from a field campaign described by Huang et al. (2011), where they collected seasonal snow in January-February 2010 on a road trip at 46 sites in North China.The sampling locations are shown in Fig. 1.At each site, samples were gathered from individual snow pits at several depths to examine snow deposited at different times during the winter.The sampling sites were chosen to be far from local sources of pollution.About 300 snow samples were collected.Processing and initial analyses were carried out at Lanzhou University in China and at three temporary laboratories, using the filtering techniques pioneered by Clarke and Noone (1985) and also used for Arctic snow by Doherty et al. (2010).Each sample was melted rapidly in a microwave oven and then immediately drawn through a 0.4 µm Nuclepore filter to extract the particulates.The BC masses in snow samples are estimated using the wavelength dependence of the measured absorption from the integrating sandwich sphere (ISSW) spectrophotometer (Grenfell et al., 2011).The spectrally-resolved total light absorption is divided into BC and non-BC fractions based on the absorption Ångström exponent of the material on the filter, and by assigning different absorption Ångström exponents to BC and non-BC light absorbing aerosols (Wang et al., 2013).The measured BC mass is an equivalent mass assuming a 550 nm BC mass absorption efficiency of 6.3 m 2 g −1 (Wang et al., 2013).This method is not possible to state with confidence that any of the light absorption in the snow is due to BC, when the soil and snow particulates in snow samples are the same color, which is particularly likely to happen when the snow is thin and patchy.For the sites (e.g., 41-46) with those snow samples, the BC mass concentration in snow is not reported (Wang et al., 2013).More details about the measurements in the North China campaign can be found in Wang et al. (2013) and Zhang et al. (2013).For comparison between simulations and observations, model results are sampled at the observation sampling time at each site in January and February of 2010.Although snow samples were gathered at several depths during the campaign, the BC mass concentrations in snow are mainly estimated at the top snow layer.Therefore, the simulated BC mass content in the top snow layer (which never exceeds 3 cm thickness) is compared with the observational values averaged at the top layer (2-5 cm depending on sites) snow samples.

CMC reanalysis
The Canadian Meteorological Center (CMC) reanalysis data set consists of Northern Hemisphere snow depth analysis data processed by the CMC.Snow depth data are obtained from surface synoptic observations, meteorological aviation reports, and special aviation reports that are acquired from the World Meteorological Organization (WMO) information system.The CMC reanalysis data set includes daily observations from 1998 at a spatial resolution of 24 km × 24 km.Monthly averages of snow depth and estimated snow water equivalent (SWE) are provided, where SWE is estimated using a density look-up table.More descriptions and the data can be found from the CMC website (http://nsidc.org/data/nsidc-0447.html).

Snowpack simulations
Before characterizing the distributions of BC and dust in snow, it is necessary to evaluate and understand potential biases of the WRF-Chem simulated snow distributions with available observations or reanalysis over North China.Figure 2 shows the spatial distributions of snow depth and SWE from the CMC reanalysis data and the WRF-Chem simulation averaged for January-February of 2010, with the North China campaign observations of snow depth at specific locations embedded in each panel.In general, the simulation and CMC reanalysis data are consistent with a spatial correlation coefficient of 0.6, both showing snow depth increases northward and with increasing topography, although the model simulates much larger snow depth near the northern boundary and snow cover extends further south.Although both model and reanalysis data capture reasonably the campaign observed spatial variation of snow depth (Fig. 2), the model overestimates the observed snow depth, while reanalysis underestimates it (Fig. 3).Overall, the reanalysis better captures the observed spatial variability with a correlation coefficient of 0.73 compared to 0.51 for the simulation, but the model matches the observed average snow depth over the 46 observation sites better, as indicated by the mean values of 13.6 and 9.7 cm from the simulation and reanalysis, respectively, compared to 13.2 cm from observations.In addition, the spatial distributions of SWE from the CMC reanalysis data and the WRF-Chem simulation averaged for January-February of 2010 are also shown in Fig. 2. Again, the model results and CMC reanalysis data show consistent spatial distribution of SWE with a correlation coefficient of 0.66.Model simulates a much larger SWE than the CMC reanalysis near the north boundary.The comparison is consistent with that for snow depth.
One key factor affecting the snow simulation is surface temperature.Figure 4   noticeable over areas north of 40 • N. Another key factor affecting the snow simulation is surface precipitation of snow.However, there is no publicly available observation for evaluation.

BC concentrations in snow
Figure 5 shows the spatial distribution of BC mass concentration in the top snow layer from the WRF-Chem simulation over North China in January-February of 2010.The campaign observed BC mass concentrations in the top snow layer (BCS) are also shown as the overlaid colored circles.The model simulates the highest BCS with > 5000 ng g −1 over northeast China (110-130 • E, 38-48 • N), where both BC deposition and snow depth are high (Figs. 2 and 6).BCS reduces away from northeast China and reaches < 50 ng g −1 near the north boundary of the domain.This is consistent with the North China campaign data that show high BCS in heavily industrialized regions and smaller concentrations (30-50 ng g −1 ) farther north in North China (51 • N). Figure 6 shows the spatial distribution of dry and wet deposition fluxes of BC on snow from the WRF-Chem simulation in January-February of 2010.The deposition fluxes are averaged spatially and temporally over only snow-covered surfaces, representing the deposition on snow.Dry and wet deposition fluxes of BC are more comparable over Central China (105-120 • E, 30-40 • N), but dry deposition dominates the total BC deposition in northeast China, where the highest BCS is found.The reason is that, over the surface that is frequently covered by snow, dry deposition is more ef- ficient than wet deposition that is limited to rainy or snow days (figures not shown).Wet deposition mainly results from resolved stratiform precipitation that accounts for more than 95 % of the total precipitation during the study period (figures not shown).
BCS are compared more explicitly with observations at each site in Fig. 7, following the sequence of each stop in the campaign from left to right on the x axis.During the campaign, the road trip started from Inner Mongolia (sites 1-15), which is relatively clean compared to northeast China.Most of the observed BCS over this region are around hundreds of ng g −1 .The cleanest snow (tens of ng g −1 ) was sampled near the north boundary of China (sites 21-24).Polluted snow was sampled in the industrial regions of northeast China (sites 26-40).There was no BCS data available for sites 41-46.The WRF-Chem simulation captures the observed sharp increase of BCS towards more polluted sites, from 30-50 ng g −1 at 51 • N to over thousands of ng g −1 at 42 • N. The model generally agrees well with the observations and has negative biases (within a factor of two) in the relatively clean sites (e.g., sites 1-32) with a median modelto-observation ratio of 1.03.It is noteworthy that the simulation exhibits large daily and diurnal variation of BCS as denoted by the box and the whisker bar, respectively, in Fig. 7.This large temporal variation of BCS is partly driven by the evolution process of snow and its BC content in the model.When snow starts accumulating on the ground, the BC mass in snow is much less than the snow mass (minimum BCS value).Then, BCS increases with the accumulation of snow due to both dry and wet deposition.When snow starts melting, BCS keeps increasing due to dry deposition till the snow disappears.Nevertheless, we would also like to point out that this large temporal variation of simulated BCS might be sometimes due to the artifact introduced by the discretization of snowpack into layers for numerical solutions.As the top snow layer becoming thinner and thinner, it retains most of its BC but losing its snow water, and therefore results in large BCS values.We therefore caution the comparison of simulated BCS with the measurements that sampled 2-5 cm snowpack as the top snow layer.
Although the simulated BCS at the sampling time at site 6 significantly overestimates the observed BCS, the simulated mean value of January-February is closer to the observation and has a negative bias.Other than site 6, differences between the simulated values at the snow sampling time and the monthly average exist in many other sites to different degrees.Due to the large temporal variability of seasonal snow and the uncertainty in simulating the exact timing of weather systems and aerosol emissions and deposition, the comparison between long-term averages of observed and simulated BCS is desired.However, the snow sampling at all sites from this North China campaign was only made at specific times relatively far apart along the road trip.In addition, more comparisons have also been conducted using the average of model results sampled within a few (e.g., ±2) hours and a few (e.g., ±2) days around the observation sampling time at each site.The comparisons do not change significantly (figures not shown).The large temporal variability of BCS may also indicate the caveats in comparing model simulated monthly mean values with observations from seasonal snow sampling at a specific time, a common practice in global modeling studies (e.g., Huang et al., 2011;Qian et al., 2014).The relative biases increase at sites 32-40 located within the industrial source regions.The model generally overestimates the BCS in these sites.It is interesting to see from Fig. 7b that the model generally has smaller biases in regions out- side of the industrial source regions but large positive biases are common near or inside the industrial source regions.It may be partly due to the sampling strategy that sites were chosen to be far from local sources (Huang et al., 2011).It may be difficult for the model to resolve the sub-grid variability, which is typically larger over the source regions with large emissions than relatively clean areas.
Model biases in simulating BCS may result from the biases in simulating BC lifecycle in the snow and atmosphere.For example, the biases may be partly due to the uncertainties in treating BC scavenging by snowmelt.As discussed above, the aerosol scavenging ratios in snow are prescribed in this study.These values are highly uncertain, and there are no direct measurements over the study region to constrain these model parameters.On the other hand, the biases may also result from the biases in simulating the atmospheric BC.Since observations of atmospheric BC concentrations are not publicly available in this region, satellite retrieved AOD is used to evaluate the general model performance in simulating aerosols in the atmosphere.The AOD retrieved from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument (Kaufman et al., 1997) on board the NASA Aqua platforms and the Multi-angle Imaging SpectroRadiometer (MISR) instrument (Diner et al., 1998(Diner et al., , 2001;;Martonchik et al., 2004) on board the NASA Terra platform are used.The MODIS "Deep Blue" product (Hsu et al., 2006) that can retrieve AOD over bright surfaces is used.MISR can also retrieve aerosol properties over highly reflective surfaces.
Figure 8 shows the spatial distribution of AOD over North China in January-February of 2010 from the Terra-MISR and Aqua-MODIS retrievals and the corresponding WRF-Chem simulation.The model results are sampled at the Terra-MISR (∼ 10:45 LT) and Aqua-MODIS (∼ 13:30 LT) overpass time, respectively, at the locations where retrievals are available.In general, MISR and MODIS show consistent spatial patterns of AOD over North China.But the magnitude of MODIS retrieved AOD is higher than that of MISR over land.Part of the difference between the MODIS and MISR retrievals may be due to their different overpass time, as aerosol concentrations may build up from morning to early afternoon.Ge et al. (2010) found that the MISR retrieved AOD has more reasonable agreements with ground-based observations over China compared to the MODIS "Deep-Blue" products.Overall, the model captures well both satellite retrievals, showing high AOD over the industrial regions of Central and northeast China.The model results are more consistent with the MISR retrievals.The region with high BC deposition fluxes (Fig. 6) also has large AOD.Although satellite AOD can generally evaluate the model performance in simulating aerosols, more explicit measurements are needed to evaluate the simulated BC in the atmosphere.In addition, some biases of simulated BCS may also come from uncertainties of BC deposition rate, particularly for dry deposition rate that dominates the BC deposition over northeast China where the simulated BCS show relatively large biases.In order to identify the uncertainty sources of BCS, concurrent measurements of atmospheric and snow BC and the deposition fluxes of BC are needed in future campaign.
Besides uncertainties in the model, the difference between observations and simulations could also result from uncertainties in the BCS estimates (Wang et al., 2013).The mea- sured BC mass is an equivalent mass with assumptions of BC optical properties in snowpack, such as assuming a 550 nm BC mass absorption efficiency of 6.3 m 2 g −1 (Wang et al., 2013).Therefore, the observations also show large uncertainty ranges (Fig. 7).In general, most model results fall into the uncertainty ranges of the observations.Furthermore, the observations at each site are assumed to be representative of the average over a 36 km × 36 km model grid cell, which may not be true in some cases due to the inhomogeneous spatial distribution of snow and BCS.

Dust concentrations in snow
Figure 9 shows the spatial distribution of dust mass concentrations in the top snow layer from the WRF-Chem simulations over North China in January-February of 2010.The dust mass concentrations in the top snow layer (DSTS) are highest over northwest China and Mongolia (90-110 • E and 40-45 • N) and Liaoning province, which generally follow the dust source distribution (Fig. 1).DSTS can reach ∼ 5 mg g −1 near the dust source regions, which is comparable to the magnitudes of the observed dust mass concentrations in snow after dust storms in the upper Colorado River basin of the USA (Painter et al., 2012).DSTS gradually reduces with distance away from the dust source regions and reaches ∼ 10 µg g −1 over Northeast and Central China and ∼ 100 ng g −1 near the northern boundary of the domain.The spatial distribution of DSTS follows the coincidence of dust deposition (Fig. 10) and snow coverage (Fig. 2).Both the MODIS and MISR retrievals show a band of moderate AOD in the northwest of the domain (Fig. 8), which is captured by the model.This band of moderate AOD is mainly contributed by dust (figures not shown) emitted locally from the Gobi Desert or transported eastward from the Taklimakan Desert to the west (Chen et al., 2014).The region with high dust deposition fluxes (Fig. 10) coincides well with this AOD band.
Although, there is no direct measurement of DSTS in the North China campaign, the simulated DSTS spatial distribution is consistent with the campaign finding that DSTS was the dominant absorbing aerosol in snow over the grasslands of northwest China and Mongolia (site 1-15), but is reduced significantly eastward as mentioned by Huang et al. (2011) and Wang et al. (2013).In addition, the North China campaign also found a significant amount of dust along with BC in snow in Liaoning province (site 40) (Huang et al., 2011;Wang et al., 2013), and DSTS is much higher than BCS near the dust source regions but comparable to, if not lower than, BCS over northeast China and Central China, where BC emissions dominate.The estimated absorbing aerosol concentrations in snowpack range from 100 to 4300 ng g −1 BC equivalent in the dust-dominated regions (e.g., Inner Mongolia), higher than those in the BC-dominated regions (e.g., northeast China) that ranged from 40 to 1600 ng g −1 (Huang et al., 2011;Wang et al., 2013).

Direct radiative forcing of BC and dust in snow
The single scattering co-albedo (1-SSA) of BC, dust, and snow grains are calculated offline following the methods described in Sect.2.2.2 and shown in Fig. 11.BC has the strongest absorption at all the five wavebands compared to dust and snow grains.BC absorption at the five wavebands decreases with increasing size from 0.05 to 0.5 µm (figures  not shown).On the contrary, absorption of dust and snow grains increases with increasing size, except for the waveband of 1.5-5.0µm, at which dust absorption is strongest at the smallest size (submicron particles).Dust and snow grains also have larger variation of co-albedo among the five wave bands than BC.Generally, absorption of snow grains increases with increasing wavelength, while that of dust decreases with increasing wavelength except for submicron dust particles.Although dust is much more absorbing than snow grains at wavelengths shorter than 1.0 µm, it is noteworthy that dust co-albedo could be smaller than that of snow grains, depending on the snow and dust effective radius and wavelengths.This indicates that dust in snow could increase snow albedo at some specific cases, unlike BC, which always reduces snow albedo.
Figure 12 shows the absorption enhancement ratio of BC (at 0.1 µm effective radius) mixed internally with snow grains, which is also calculated offline (Sect.2.2.2) at five spectral wavebands as a function of snow grain effective radius.BC absorption is enhanced when BC is mixed internally with snow grains compared to mixed externally, as indicated by the enhancement ratios larger than one, except for the waveband of 1.5-5.0µm, at which BC absorption could be reduced when mixed internally with snow grains larger than ∼ 250 µm effective radius.In general, the BC absorption enhancement effect increases with decreasing wavelength and decreases with increasing snow grain effective radius.In most cases, BC absorption could be enhanced by a factor of 1.5-2 dependent on the snow grain effective radius and wavelength.
Figure 13 shows the spatial distributions of direct radiative forcing by BC and dust in the atmosphere and in snowpack from the WRF-Chem simulation over North China in January-February of 2010.The forcing is averaged when snow is present, representing the mean increase in energy absorption by snowpack.It is shown that both BC and dust introduce radiative warming in snow.The warming patterns are consistent with the spatial distributions of BCS and DSTS (Figs. 5 and 9).The warming effect reaches ∼ 10 W m −2 for both BC and dust in snow.Although dust is much less absorbing than BC (Fig. 11), it results in similar magnitude of warming in snow as BC due to its much larger mass content in snow than BC.As a reference, the direct radiative forcing of atmospheric BC and dust at the surface is also shown.The patterns of direct radiative forcing of atmospheric BC and dust follow the distributions of their concentrations in the atmosphere (figures not shown).Atmospheric BC and dust reduce radiation reaching the surface and have a cooling effect of ∼ −10 W m −2 near the source regions.The pattern of BC warming in snow shows higher magnitude in the south than in the north, while that of BC cooling at the surface is reversed.The cooling pattern at the surface of atmospheric dust roughly coincides with the warming pattern induced by dust in snow.It is interesting to note that BC and dust in snow result in comparable magnitudes of surface radiative forcing as in the atmosphere but with an opposite sign, which may indicate that BC and dust radiative forcing in snow and at the surface could offset each other in winter over some regions of North China, particularly for dust.

Conclusions
In this study, the WRF-Chem model is coupled with the SNICAR model to simulate the BC and dust concentrations and their radiative forcing in seasonal snow over North China in January-February of 2010, consistent with a field campaign region and period.This study represents the first effort in evaluating WRF-Chem for simulating seasonal snowpack and its impurities against various observations made over North China.The direct radiative forcing of BC and dust in snow is also estimated for the first time at a relatively high spatial resolution over North China.In general, the model captures well the observed spatial distributions of surface temperature, snow properties, and aerosol contents in snow.The difference between observed and simulated BCS may be partly because the observations at each site are assumed to be representative of the 36 km × 36 km model grid cell, which may not be true in some cases due to the inhomogeneous spatial distribution of snow and BCS.The simulated large daily and diurnal variation of BCS cautions comparison of model simulated mean values with observations from seasonal snow sampling at a specific time adopted in some global modeling studies (e.g., Huang et al., 2011;Qian et al., 2014).
The model simulates similar magnitudes of BCS and DSTS induced radiative warming in snow.The magnitudes of radiative warming of BCS and DSTS are comparable to those of the surface radiative cooling due to BC and dust in the atmosphere, which may indicate that BC and dust radiative forcing in snow and at the surface could offset each other in winter over some regions of North China.Our estimated BC forcing (∼ 10 W m −2 ) in snow is larger than the annual mean value (∼ 1-4 W m −2 ) (averaged spatially and temporally over only snow-covered surface) estimated by Flanner et al. (2007) for North China.One reason could be that our estimation is for 2010 that has more than double the BC emissions in 1998 used in Flanner et al. (2007).Another reason could be the difference of period for average (two months versus the snow seasons of entire year).It may also be partly due to model differences.Our estimated dust radiative forcing in snow is roughly consistent with that in the upper Colorado River basin (tens of W m −2 ) induced by several mg g −1 dust in snow (Skiles et al., 2012).Our estimated dust radiative forcing per unit dust mass in snow is smaller than that of Skiles et al. (2012), perhaps partly due to our estimation being for January-February, when the surface insolation is weaker than that of Skiles et al. (2012) for March-April.
There are also many uncertain factors and processes affecting the simulation of aerosols and their radiative forcing in snowpack.Besides aerosol refractive index, shape, and mixing state with other materials discussed in Flanner et al. (2012), aerosol size distribution is also critical, not only in snow but also in the atmosphere.The size distribution is particularly important for dust, since it could significantly affect the dust deposition fluxes (Zhao et al., 2013c) and determine the relative absorption of dust with respect to snow grains.Unlike BC, dust could be less absorbing than snow grains at near-infrared wavebands depending on the sizes.In one of our sensitivity simulations with the modal approach to represent dust size distributions as in Zhao et al. (2013c), dust could result in cooling effect in snow in some specific cases.In addition, following Flanner et al. (2012), the size-resolved atmospheric BC is lumped together in snow and assumed a fixed effective radius.More work is needed to resolve the BC size distribution in snow.Furthermore, the atmospheric dust, if mixed internally with other aerosols, can be activated as cloud droplets.However, the effect of particles mixed internally with snow grains is only considered for BC in snowpack.Dust particles are all assumed to mix externally with snow grains.BC amounts in snow grains may also be underestimated.Although this study accounts for the ice nucleation of liquid droplets, which may contain BC through cloud droplet activation, ice nucleation of aerosols for iceand mixed-phase clouds are not included.Therefore, the deposition of ice-borne aerosols is not considered.Finally, the inconsistency between the size distributions of deposited dust and that used for dust optical properties calculation in snowpack may also introduce uncertainties.Future study should define consistent optical properties for the atmospheric and snow dust.
In this study, although the simulated spatial distributions of BC and dust content in snowpack are generally evaluated with observations, it is still difficult to fully quantify and attribute the model biases in simulating snow impurities, which could result from uncertainties in many model processes, including emissions, deposition processes, and snow melting scavenging.To appropriately represent BC and dust content in snow, the model needs to simulate reasonably well the aerosol lifecycle in snow as discussed by Flanner et al. (2007Flanner et al. ( , 2012) ) and Qian et al. (2014).The parametric uncertainties of the SNICAR model need to be quantified and constrained with observations.Furthermore, the aerosol lifecycle in the atmosphere should also be reasonably well simulated.The model needs to reproduce the atmospheric aerosol concentrations and treat the deposition processes correctly.Although the WRF-Chem simulated AOD is generally consistent with satellite retrievals, no direct measurement of BC or dust in the atmosphere is available for comparison.In addition, there is also no measurement to evaluate the simulated deposition fluxes.In order to fully understand model biases in simulating snow impurities, concurrent measurements of atmospheric and snow aerosols and the deposition fluxes of aerosols are needed in future field campaigns.
Finally, this study only estimates the direct radiative forcing by BC and dust in snowpack, which results from the enhanced absorption of solar radiation by BC and dust.It should be noted that BC and dust can indirectly enhance snow absorption by increasing the grain size due to acceleration of grain growth from the direct effect (first indirect effect) and by the earlier exposure of a darker substrate (second indirect effect) (Flanner et al., 2007), which cannot be diagnosed in a single experiment.Further studies are needed to examine BC and dust impact with both direct and indirect effects in snow and hence the impact on climate and hydrological cycle, which can be quantified with numerical experiments by including and excluding aerosols in snow.The WRF-Chem model coupled with SNICAR can also be applied to other regions of the globe to understand the impact of BC and dust in snowpack on the regional climate and water supplies.

Figure 1 .
Figure 1.Model domain with BC and dust emissions in January-February 2010.The 46 observation sites in the North China campaign are denoted by pentagrams in black.
at 0.5 • × 0.5 • horizontal resolution for 2006 except that BC, OM, and sulfate emissions over China are from the China emission inventory for 2010 described by Lu et al. (2011) at a 0.1 • × 0.1 • horizontal spatial resolution and a monthly temporal resolution for the simulation period.Biomass burning emissions for the simulation period (October 2009-February 2010) are obtained from the Global Fire Emissions Database, Version 3 (GFEDv3) with a monthly temporal resolution and a 0.5 • × 0.5 • horizontal resolution (van der Werf et al., 2010) and vertically distributed following the injection heights suggested by Dentener et al. (2006) for the Aerosol InterComparison project (AeroCom).Sea salt and dust emissions follow Zhao et al. (2013c).Vertical dust fluxes are calculated with the the Georgia Institute of Technology-Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GO-CART) dust emission scheme

Figure 2 .
Figure 2. Spatial distribution of snow depth (cm) and snow water equivalent (cm) from the WRF-Chem simulations and CMC reanalysis over North China in January-February 2010.Snow depth (cm) from the campaign observations (colored circle) is also shown.
shows the spatial distribution of 2 m temperature from the observations (colored circle) and the corresponding WRF-Chem simulation over North China in January-February of 2010.Daily observations of 2 m temperature in January-February of 2010 at 225 sites in North China are obtained from the Chinese National Meteorological Center (CNMC).The simulation can capture well the observed 2 m temperature with spatial correlation coefficient of 0.95 over the region of interest.The 2 m temperature decreases dramatically from around freezing level to about −30 • C with increasing latitude.The simulated 2 m temperature averaged over the observation sites is −13 • C, which is 2 • C colder than that from observations, but the caveats in comparing grid mean simulated values with station data should be noted.The model cold bias is particularly

Figure 3 .
Figure 3.The campaign observations of snow depth (cm) versus the corresponding WRF-Chem simulations (black dots) and CMC reanalysis (red dots).

Figure 4 .
Figure 4. Spatial distribution of 2 m temperature from the observations (colored circle) and the corresponding WRF-Chem simulations over North China in January-February 2010.

Figure 5 .
Figure 5. Spatial distribution of BC mass concentration in the top snow layer from the WRF-Chem simulations over North China in January-February 2010.The campaign observed BC mass concentrations in the top snow layer (colored circle) are also shown.

Figure 6 .
Figure 6.Spatial distribution of dry and wet deposition fluxes of BC on snow from the WRF-Chem simulations over North China in January-February 2010.

Figure 7 .
Figure 7. (a) BC mass concentrations in the top snow layer from the campaign observations (black) and the corresponding WRF-Chem simulations (red).The box plot of observations shows the minimum and maximum possible values from estimates and the bar within the box shows the most likely value; the whisker bar of model results shows the 10th and 90th percentiles of simulated hourly values of January-February.The asterisk represents the 24 h average value for January-February.The box of model results shows the 10th and 90th percentiles of simulated values at observation sampling hour of January-February.The bar within the box represents the simulated value at the observation sampling hour and day.(b) Ratio (red circles) of BC mass concentration in the top snow layer from the model at the observation sampling time and observation along with latitudes (green line).

Figure 8 .
Figure 8. Spatial distribution of AOD over North China in January-February 2010 from the MISR and Aqua-MODIS retrievals and the corresponding WRF-Chem simulations.The model results are sampled at the MISR (∼ 10:45 LT) and Aqua-MODIS (∼ 13:30 LT) overpass time, respectively, at the locations where retrievals are available.

Figure 9 .
Figure 9. Spatial distribution of dust mass concentration in the top snow layer from the WRF-Chem simulations over North China in January-February 2010.

Figure 10 .
Figure 10.Spatial distribution of dry and wet deposition fluxes of dust on snow from the WRF-Chem simulations over North China in January-February 2010.

Figure 11 .
Figure 11.Single scattering co-albedo (1-SSA) at the five spectral wavebands of snow grains as a function of effective radius (colored solid lines), BC with an effective radius of 0.1 µm (asterisk), and dust with four different sizes (Dust1: circle, Dust2: square, Dust3: triangle, and Dust4: diamond).

Figure 12 .
Figure 12.Absorption enhancement ratio of BC (0.1 µm effective radius) internally mixed with snow grains at the five spectral wavebands as a function of snow grain effective radius.

Figure 13 .
Figure 13.Spatial distribution of BC and dust direct radiative forcing at the surface and in the snow from the WRF-Chem simulations over North China in January-February 2010.

Table 1 .
Matching of dust size distributions between WRF-Chem and SNICAR.The values represent mass fractions of deposited atmospheric dust from WRF-Chem into each size of dust in snow in SNICAR.