Interactive comment on “Impacts of absorbing aerosol deposition on snowpack and hydrologic cycle in the Rocky Mountain region based on variable-resolution CESM (VR-CESM) simulations”

This is a very end-to-end modeling analysis of the effects on surface air temperature and snowpack (SWE, fraction and runoff) in several regions of the western U.S. It includes a comparison of modeled near-surface atmospheric BC concentrations and mixing ratios of BC and dust in snow against observational datasets. I have no fundamental problems with the analysis. The paper should be accepted after addressing the issues noted below. It could use some editing for English but overall is well-written, if a bit long, in part due to being repetitive in some places. I have enclosed an annotated version of the paper showing the small edits I think are needed for better


Introduction
Water resources are essential to human society and economic development as well as ecosystems in the western United States.Most primary water resources in the inland western US come from the Rocky Mountain snowpack (Serreze et al., 1999).Therefore, to develop a water resource management strategy, it is necessary to have information on snow accumulation and snowmelt timing.Climate change is an important factor influencing the snowpack in the Rocky Mountain region, as has been shown in many previous studies (e.g., Abatzoglou, 2011;Pederson et al., 2011;Rhoades et al., 2017).Another important factor is the light-absorbing aerosols (LAAs, e.g., black carbon (BC), organic carbon (OC), and dust) in snow (e.g., Flanner et al., 2007;Painter Published by Copernicus Publications on behalf of the European Geosciences Union.The refined region at a resolution of 0.125 • is surrounded by dashed lines.(c) Five regions are identified for the analysis in this study, including three mountainous region (1, Northern Rockies; 2, Greater Yellowstone region; 3, Southern Rockies) and two regions in the plains around the mountains (4, eastern Snake River Plain; 5, southwestern Wyoming).et al., 2007;Qian et al., 2015;Yasunari et al., 2015).Previous studies have shown that LAAs in snow can significantly reduce the surface albedo (often known as the snow darkening effect; SDE), modify the surface energy budget and snowmelt, and lead to the modification of hydrologic cycles (e.g., Warren and Wiscombe, 1980;Hansen and Nazarenko, 2004;Flanner et al., 2007Flanner et al., , 2009;;Painter et al., 2007Painter et al., , 2010;;Qian et al., 2009Qian et al., , 2011;;Yasunari et al., 2015).Moreover, the LAA-induced snow albedo reduction may initiate positive feedback processes, which can amplify the reduction of snowpack (e.g., Flanner et al., 2009;Qian et al., 2009).
In past decades modeling studies have been undertaken to quantify the impacts of SDE by LAAs (e.g., Flanner et al., 2007;Qian et al., 2009;Oaida et al., 2015;Yasunari et al., 2015).Generally the models developed have the ability to simulate the temporal evolution of snow albedo under the influence of LAAs in snow.These studies have enhanced our understanding of the spatial and temporal variations in climate forcings and responses due to LAAs in snow from regional scales (e.g., Qian et al., 2009;Oaida et al., 2015) to global scales (e.g., Flanner et al., 2007;Yasunari et al., 2015).
For example, the impacts of LAAs in snow are stronger in regions with considerable snow cover and sufficient LAA deposition (e.g., Arctic, Northeast China, Tibetan Plateau, and western US), and they are largest during the snowmelt period due to the positive snow albedo feedback.However, as also mentioned in these studies, a reliable quantification of the impacts of LAAs in snow is hindered by the model deficiencies in simulating the snowpack and aerosol cycles, with additional uncertainties induced by the parameterization of snow-aerosol-radiation interactions.
In particular, previous studies have used coarse-resolution global climate models (GCMs) or high-resolution regional climate models (RCMs) to quantify the impacts of LAAs in snow.However, there are weaknesses both in coarseresolution GCMs and in RCMs.Both snowfall and snow accumulation depend on temperature and precipitation, and thus the distribution of snowpack depends strongly on topographic variability.Current GCMs with a typical horizontal resolution of 1 to 2 • cannot resolve the snowpack over regions with complex terrains (e.g., Rocky Mountains) due to the coarse resolution (Rhoades et al., 2016;Wu et al., 2017), which impedes the reliable quantification of SDE by LAAs in mountainous regions (e.g., Flanner et al., 2007;Yasunari et al., 2015).RCMs can simulate the snowpack more accurately than coarse-resolution GCMs, but they are not able to simulate the global transport of aerosols to the focused region except when aerosol transport along the boundary is prescribed (e.g., Qian et al., 2009).Moreover, LAAs in snow may also influence the climate beyond the focused region (e.g., Yasunari et al., 2015), which cannot be accounted for in RCMs.Variable-resolution GCMs (VR-GCMs) can overcome these weaknesses of either coarse-resolution GCMs or RCMs and serve as a better tool to quantify the impacts of LAAs in snow.Although GCMs with globally uniform high resolutions (10-30 km) may be an ideal tool to simulate the snowpack and snow-aerosol-radiation interactions, they are not widely applied due to the constraints of computational resources (e.g., Haarsma et al., 2016).Instead, using VR-GCMs is a more economic approach and has gained increasing utility in recent years (e.g., Zarzycki et al, 2014a, b;Sakaguchi et al., 2015).

A variable-resolution version of the Community Earth
System Model (VR-CESM) has been developed (Zarzycki et al., 2014a, b).With a refined high resolution, VR-CESM has shown significant improvements of the Atlantic tropical storms (Zarzycki and Jablonowski, 2014) and South America orographic precipitation (Zarzycki et al., 2015).The model has also been used in the regional climate simulations over the western US, and the results show that the VR-CESM is capable of reproducing the spatial patterns and the seasonal evolution of temperature, precipitation, and snowpack in the Sierra Nevada (Huang et al., 2016;Rhoades et al., 2016) and Rocky Mountains (Wu et al., 2017).In particular, VR-CESM reasonably simulates the magnitude of snow water equivalent, the timing of snow water equivalent peaks, and the duration of snow cover in the Rocky Mountains, as shown in comparison with Snow Telemetry (SNOTEL) and MODIS (Moderate Resolution Imaging Spectroradiometer) snow cover observations (Wu et al., 2017).
Following the evaluation study of Wu et al. (2017), here we use VR-CESM to investigate the impacts of LAAs in snow (BC and dust) on the snowpack and hydrologic cy-   (Skiles and Painter, 2016b) cles over the Rocky Mountains.By comparing the two VR-CESM simulations with and without LAAs in snow, we examine the impacts on surface radiative transfer, temperature, snowpack, and runoff induced by LAAs in snow.To our knowledge, it is the first time that VR-CESM is applied for the study of LAAs in snow.Our results will demonstrate that VR-CESM is skillful for this kind of research.The remainder of the paper is organized as follows.Section 2 introduces the model and experimental design.Section 3 describes the observation data used for the validation of model simulations of aerosol fields in the surface air and in snow.Section 4 presents the evaluation of aerosols fields, followed by their surface radiative effect (SRE), and the change in surface temperature, snowpack, and runoff induced by LAAs in snow.Discussion and conclusions are given in Sect. 5.

Model and experimental design
The model used in this study is VR-CESM, a version of CESM (version 1.2.0) with variable-resolution capability (Zarzycki et al., 2014a, b).CESM is a state-of-the-art Earth system modeling framework that allows for investigation of a diverse set of Earth system interactions across multiple timescales and space scales (Hurrell et al., 2013).CESM uses the Community Atmosphere Model version 5 (CAM5) for the atmospheric component (Neale et al., 2010).The variable-resolution capability is implemented into the spectral element (SE) dynamic core of CAM5.The SE dynamic core uses a continuous Galerkin spectral finiteelement method designed for fully unstructured quadrilateral meshes and has demonstrated near-optimal (close to linear) parallel scalability on tens of thousands of cores (Dennis et al., 2012).This enables the model to run efficiently on decadal to multi-decadal timescales.For the land component, CESM uses the Community Land Model version 4 (CLM4).CLM4 can be run at the same horizontal resolutions as CAM5 and can thus also benefit from the variableresolution capability of CAM5.
CESM also includes advanced physics for CAM5 (Neale et al., 2010) and CLM4 (Oleson et al., 2010).The CAM5 physics suite consists of shallow convection (Park and Bretherton, 2009), deep convection (Zhang and McFarlane, 1995;Richter and Rasch, 2008), cloud microphysics (Morrison and Gettelman, 2008) and macrophysics (Park et al., 2014), radiation (Iacono et al., 2008), and aerosols (Liu et al., 2012).For aerosols, a modal aerosol module (MAM) is adopted to represent the internal and external mixing of aerosol components such as BC, OC, sulfate, ammonium, sea salt, and mineral dust (Liu et al., 2012).Here, we use the three-mode version of MAM (MAM3).These three modes are Aitken, accumulation, and coarse modes.In MAM3, BC is treated in the accumulation mode.BC particles are instan-taneously mixed with sulfate and other components in the accumulation mode once emitted.Dust particles with the diameter range of 0.1-1 and 1-10 µm are emitted into the accumulation mode and coarse mode, respectively.Airborne aerosol particles are then transported by winds and delivered back to the land surface by both dry and wet deposition, as described in Liu et al. (2012).
CLM4 physics includes a suite of parameterizations for the land-atmosphere exchange of water, energy, and chemical compounds.In particular, CLM4 explicitly represents the snowpack (accumulation due to snowfall and frost, loss due to sublimation, and melt) by a snow model and its coupling with the SNow, ICe and Aerosol Radiation (SNICAR) model for snow-aerosol-climate interactions (Flanner et al., 2007).SNICAR incorporates a two-stream radiative transfer solution of Toon et al. (1989) to calculate the snow albedo and the vertical absorption profile from the solar zenith angle, the albedo of the substrate underlying snow, mass concentrations of atmospheric-deposited aerosols (BC and dust), and ice effective grain size (r e ); r e is simulated with a snow aging routine (Oleson et al., 2010).SNICAR is compatible with the new modal aerosol module of CAM5 in the treatment of aerosol deposition (Liu et al., 2012).It should be mentioned that SNICAR includes the effects of feedbacks to the snowpack (grain size, melt) that are driven by the snow albedo reduction due to LAA deposition.As our knowledge of OC optical properties is limited, the impact of absorbing OC on snow albedo is not included in the standard CLM4 and thus not considered in this study.Note that the SNICAR model we use assumes spherical snow grains and aerosolsnow external mixing for the calculation of snowpack optical properties (Flanner et al., 2007;Oleson et al., 2010).Recent studies have shown that non-spherical snow grains play a critical role in snow albedo calculations and decrease the snow albedo reductions included by LAAs compared with spherical snow grains (e.g., Liou et al., 2014;Dang et al., 2016;He et al., 2014He et al., , 2017)).Nonetheless, the knowledge of snow grain shape evolution is limited and thus spherical snow grains are assumed.Studies have also shown the significant enhancement of solar radiation absorption with larger snow albedo reductions by aerosol-snow internal mixing compared to aerosol-snow external mixing (e.g., Flanner et al., 2012;He et al., 2014;Liou et al., 2014).However, without considering aerosol-snow internal mixing, the SNICAR model we use assumes absorption-enhancing sulfate coatings to hydrophilic BC, which can mimic BC coatings by snow and compensate for the neglect of the absorption enhancement by aerosol-snow internal mixing (Flanner et al., 2007(Flanner et al., , 2012)).Therefore, the impacts of BC-in-snow shown in this study (Sect.4) are not necessarily biased low.Despite this, by assuming dust-snow external mixing this study may underestimate the impacts of dust-in-snow.
For the high-resolution modeling, we have designed a variable-resolution grid that transits from global quasiuniform 1 • resolution to a refined 0.125 • resolution in the Rocky Mountains (Fig. 1a).The variable-resolution grid is the same as that used in Wu et al. (2017) and is generated by the open-source software package called SQuadGen (Ullrich, 2014).A topographical dataset for this variable-resolution grid is also generated accordingly by the National Center for Atmospheric Research (NCAR) global model topography generation software called NCAR_Topo (v1.0) (Lauritzen et al., 2015) as described in Wu et al. (2017).Figure 1b shows the spatial variations in terrain height for the variable-resolution grid used in VR-CESM.Compared to United States Geological Survey (USGS) 3 km topography data (Lauritzen et al., 2015), the topography data used in VR-CESM resolve the variations in terrain in the Rocky Mountains well (see Fig. 2 of Wu et al., 2017).In Wu et al. (2017), we have shown that VR-CESM performs well in the simulation of regional climate patterns, including spatial distributions and the seasonal evolution of temperature, precipitation, and snowpack in the Rocky Mountain region.In this study, we further apply VR-CESM to simulate the SDE of LAAs and its impacts on snowpack and hydrologic cycles in the Rocky Mountains.
VR-CESM is run in the coupled land-atmosphere mode with prescribed observed monthly 1 • × 1 • sea surface temperature and sea ice coverage (Hurrell et al., 2008), following the Atmospheric Model Intercomparison Project (AMIP) protocols (Gates, 1992).The simulation period is from 1979 to 2005, and the results for the last 25 years  are used for the analysis shown below.Historical greenhouse gas concentrations and anthropogenic aerosol and precursor gas emissions are prescribed from the datasets of Lamarque et al. (2010).In particular, the BC emissions consist of various sources, including domestic, energy, transportation, waste, shipping, and wildfire (forest and grass fires) emissions.The horizontal resolution for BC emission used in this study is 1.9 × 2.5 • .We note that BC emission data are natively at a resolution of 0.5 • × 0.5 • (Lamarque et al, 2010).However, they are processed to be at a relatively coarse resolution of 1.9 • × 2.5 • for adoption in standard CESM, which is used in this study.The relatively coarse resolution of BC emission may partly explain the model bias in the simulation of BC concentrations near the surface and in snow across regions where local BC sources can contribute significantly to the observed BC concentrations, as will be discussed in Sect. 4. It is desirable to adopt BC emission at its native resolution for our high-resolution simulation.The sensitivity of our simulation results to the resolution of BC emission will be analyzed in a separate study.
For dust aerosol, the emission flux is calculated interactively in the model at each time step by using a dust emission scheme (Oleson et al., 2010).The dust emission flux is calculated from the friction velocity, threshold friction velocity, atmospheric density, clay content in the soil, areal fraction of exposed bare soil, and source erodibility (Oleson et al., 2010;Wu et al., 2016).Due to the large uncertainty in modeled dust emission, the dust emission scheme also adopts a tuning fac- VR-CESM simulation and IMPROVE observations.Also given are the mean results at all the stations from simulation and observations and their correlation coefficient (R).The 1 : 1 (solid) and 1 : 5/5 : 1 (dashed) lines are plotted for reference.
tor (T ) to simulate a reasonable dust emission amount.Our test simulation shows that with the increase in model resolution, VR-CESM produces much higher dust concentrations compared to the observations (Sect.3) in North America if the T used in the standard CESM with a quasi-uniform 1 • resolution is used.Therefore, for VR-CESM simulation in this study, T is reduced by a factor of 2.6 to produce similar magnitudes of near-surface dust concentrations as the observations, as will be shown in Sect.4.1.Note that such a reduction of T is only applied in North America, since other continents have a resolution of quasi-uniform 1 • , the same as in the standard CESM.
In addition to the control experiment with the impacts of LAAs (BC and dust) in snow included (CTL), we conduct a sensitivity experiment that turns off the impact of LAAs in snow (NoSDE).Through the comparison of these two simulations (CTL and NoSDE), the impacts of SDE by LAAs on the snowpack and hydrologic cycles can be identified.To facilitate the analysis of SDE, we also calculate the surface radiative effect (SRE) by BC-in-snow in the control experiment from the difference in absorbed radiation with all aerosols (i.e., the standard radiation call) and with all aerosols except BC (i.e., only dust in this case; a diagnostic radiation call) as in Flanner et al. (2007).SRE by dust and by BC and dust are calculated similarly.
To quantify the impacts of LAAs in snow, we mainly focus on five regions.Three of these regions are in the high mountains: Northern Rockies, Greater Yellowstone region, and Southern Rockies.The elevation is higher in the Greater Yellowstone region and Southern Rockies (> 2250 m) than in the Northern Rockies (> 750 m).The other two regions are over the plains near the mountains: Snake River Plain and southwestern Wyoming.These two regions are selected be-cause they are close to the source regions of BC and dust and also have considerable snow cover (> 50 %) in winter.These five regions are shown in Fig. 1c.

Observations
We will use various observations to validate the model simulation of aerosol (BC and dust) concentrations near the surface and in snow.
First, we use the observations of near-surface atmospheric BC and dust concentrations from the Interagency Monitoring of PROtected Visual Environments (IMPROVE) network (Malm et al., 1994).Observed mass concentrations of elemental carbon (EC) are used for comparison with model simulation of BC concentrations.Although EC can be somewhat different from BC (Andreae and Gelencser, 2006), EC concentrations have been widely used for the validation of modeled BC concentrations in previous studies (e.g., Koch et al., 2009;Liu et al., 2012).For dust, simulated dust concentration accounts for dust particles with diameters below 10 µm.To compare, observed mass concentrations of fine soil (FS, with a diameter < 2.5 µm) and coarse mass (CM, with a diameter between 2.5 and 10 µm) from IMPROVE are combined, following the approach of Kavouras et al. (2007) and Wells et al. (2007).In reality, in addition to dust, CM may also contain other aerosols such as sulfate, nitrate, organic and elemental carbon, and sea salt.However, according to the study of Malm et al. (2007), who analyzed the speciation of coarse particles collected at nine selected rural IMPROVE stations in 2004, the contributions of dust to CM are above 70 % (74-90 %) at the three stations in the inland western US.In their study, lower contributions of dust to CM (34 and 65 %) were Note that the BC and dust mass mixing ratios are given in different units, i.e., ng g −1 and µg g −1 , respectively.found in the two stations near the coast.We caution that these two stations were < 150 km away from the metropolitan regions, indicating that urban emissions may also contribute to CM there.Additional contributions may result from sea salt or sodium nitrate resulting from reactions of nitric acid with sea salt, as mentioned in their study (Malm et al., 2007).Therefore, to minimize the contributions of other aerosols to CM, we do not use the stations in or near the metropolitan regions or near the coast for the validation of dust concentration.Nonetheless, we acknowledge that there may be small contributions from other aerosols to CM, and the estimated dust concentration by summing FS and CM may represent an upper limit of dust concentrations (with a diameter < 10 µm) from the observations.Note that the observation period of IMPROVE varies with the stations; some stations started collecting data in the 1980s and some more recently (2000s).To derive a climatological dataset for model comparisons, we only select the stations with more than 5 years of dust observations.In total 80 and 94 stations are selected for BC and dust observations, respectively, in the western US (Fig. 2).
Second, we use field measurements of BC mass mixing ratio in snow (C BC ) from previously published studies.Although field observations of C BC in snow extend back to the 1980s, they were made mostly in the polar regions, the Alps, the Cascade Mountains, eastern Canada, and western Texas and New Mexico (see Qian et al., 2015 and references therein).Recently, Doherty et al. (2014) made valuable measurements of the vertical profiles of LAAs in seasonal snow from January to March 2013 in the western US.They used an integrating sphere-integrating sandwich (ISSW) spectrophotometer to estimate the C BC over 67 sites in North America (including 17 sites in the Rocky Mountain region).Observed C BC by Doherty et al. (2014) was recorded on a single day.Doherty et al. (2016) further provided the temporal The locations and sample dates as well as the measurements for these stations are given in Table 1.For comparison with model simulations, we derive observed BC mass mixing ratios (C BC ) in the whole snow column at sites 1-12 and 16-17 by dividing total BC mass throughout the snow column by total snow mass throughout the snow column.At sites 13-15, the averages of C BC for all the aged snow samples (from various depths and columns) were reported by Doherty et al. (2016) and are used here.If measurements of C BC on multiple days were made, the means and standard deviations of C BC are given.As our simulation period  does not encompass the years 2013 and 2014, we will use the daily simulation results of C BC on the same month and day (or months and days; Table 1) when the observations were made (i.e., we will ignore the exact year) and compare them (means and standard deviations) with the observations.At each station, daily simulation results are used only when snow is present (i.e., daily mean snow water equivalent ≥ 1 mm); 1 mm is chosen to be consistent with the minimum snow-layer thickness in observations.
There are few observations of dust mass mixing ratio in snow (C dust ) in the Rocky Mountain region.To our best knowledge, the only published observations were conducted at two sites in the Southern Rockies: one in the Senator Beck Basin Study Area (SBBSA) in the San Juan Mountains with at least 9-year (2005-2013) records (Painter et al., 2012;Skiles et al., 2015;Skiles and Painter, 2016a, b) and the other in the Grand Mesa (∼ 150 km to the north of SBBSA) with at least 4-year (2010-2013) records (Skiles et al., 2015).Snow samples in the top 30 cm of the snow column were collected at irregular time intervals from March to June.Here we will use the end-of-year (EOY) C dust , which was reported from the samples collected just prior to snow  1) from VR-CESM simulations and observations, with the error bars denoting the corresponding standard deviations.The observations are compiled from previously published studies (Table 1).If multiple observations are recorded at a certain site, the observed standard deviations are calculated from these multiple observations (Sect.3).Simulated BC mass concentration in the snow column and its standard deviation are calculated from the 25-year mean and standard deviation of simulation on the same month and day as the observations (Sect.3).The 1 : 1 (solid) and 1 : 5/5 : 1 (dash) lines are plotted for reference.depletion and consisted of the majority of dust in the snow column (Skiles et al., 2015;Skiles and Painter, 2016a).For the simulation, we will calculate mean C dust for May-June from daily C dust on the days when snow is present (i.e., snow water equivalent ≥ 10 mm).Another consideration is that observed C dust contains all the dust particles, while simulated C dust only accounts for the dust particles with diameters smaller than 10 µm.According to the observations by Reynolds et al. (2016), the mass concentration of total suspended particles (TSP) both in the atmosphere and in snow is mainly from particles with diameters larger than 10 µm in the Utah-Colorado region.This will affect the model comparison with the observations, which will be discussed in Sect. 4.

Spatial patterns of near-surface aerosol concentrations
Before we examine the impacts of aerosol deposition onto snow, we will first evaluate the aerosol simulations by the model.Figure 2 shows the spatial patterns of cold season (winter and spring) mean emission fluxes and near-surface concentrations of BC and dust in the western US from the VR-CESM simulation.The IMPROVE stations are also de-noted by circles, with larger circle sizes indicating higher observed near-surface BC or dust concentrations.In the model, the BC emission flux is prescribed and is largest on the Pacific Coast and in southern Arizona.BC emission fluxes are relatively large in central-northern Colorado and northwestern Utah, where large metropolises are located.Corresponding to the patterns of BC emission flux, simulated nearsurface BC concentrations are also higher (> 100 ng m −3 ) in these regions.A band with relatively high near-surface BC concentrations around 50-100 ng m −3 is also found in southern Idaho to the west of the Greater Yellowstone region and to the south of the Northern Rockies, indicating the transportation of BC around the mountains.Near-surface BC concentrations decrease at higher elevations.The spatial patterns simulated by the model are generally consistent with observations; higher BC concentrations in the source regions and lower in the mountains.
Dust sources are located in the dry regions with exposed bare soils, such as the southwestern US (southern California, western Arizona, and southern New Mexico), northern Mexico, the Great Basin, and the Colorado Plateau.Dust emissions are also found in the Great Plains, although they are much weaker.In the Great Plains, agricultural activities can disturb the soil, making it vulnerable to wind erosion (Ginoux et al., 2012).Simulated cold season mean dust concentrations are higher (10-500 µg m −3 ) in the source regions, but decrease dramatically (0.1-5 µg m −3 ) toward the mountains.Compared to the observations, the model reproduces the spatial patterns of near-surface dust concentrations with higher concentrations in the southwestern part of the US.However, the model tends to overestimate the dust concentrations in Utah, indicating that dust emissions may be overestimated there.
Comparisons of modeled and observed near-surface BC and dust concentrations at the IMPROVE stations are further shown in Fig. 3.The modeled concentrations are generally within a factor of 5 of the observed concentrations, and the two are moderately correlated (the correlation coefficients (R) are 0.56 and 0.47 for BC and dust concentrations, respectively).Averaged across all comparison stations, the modeled BC concentration is a factor of 1.8 lower than the observed concentrations, and the modeled dust concentration a factor of 1.4 higher.The model tends to systematically underestimate observed near-surface BC concentrations in the Utah-Nevada region, the Rocky Mountains, and the Great Plains, where the stations are located downwind of source regions (Fig. 2b).In particular, observed near-surface BC concentrations are underestimated mostly by a factor of 1.5-5 in the Rocky Mountains.The underestimation of nearsurface BC concentrations in these regions may suggest that the transport of BC in our simulations is too weak.This deficiency may also be ascribed to local BC sources (e.g., Doherty et al., 2014) not resolved by the prescribed BC emission in the model (e.g., at 1.9 × 2.5 • resolution).For dust, although the model overestimates near-surface dust concen- trations for most of the stations near the dust sources (southwestern US, Utah, and Nevada), the model reasonably simulates the magnitude of near-surface dust concentrations in the Rocky Mountains.This may also be associated with underestimated transport in the model, which is consistent with the low bias in near-surface BC concentrations in downwind regions.
Note that although only the BC and dust emission fluxes over the western US are shown in Fig. 2, the long-range transport of these aerosols from other regions (e.g., Asia and Africa) can also contribute to BC (e.g., Zhang et al., 2015) and dust (Wells et al., 2007) concentrations in the western US.In addition, there are substantial variations in aerosol emissions in the western US.As mentioned in Sect.2, although we adopt VR-CESM with a refined high resolution (0.125 • ) in the Rocky Mountains, we use a coarse-resolution gridded emission dataset (i.e., 1.9 • × 2.5 • ) for BC.For dust, the small-scale variations in dust emissions can be represented in the model as they are calculated online in the model.However, dust emission depends on many variables such as near-surface winds, soil moisture, vegetation cover, and soil texture (Oleson et al., 2010;Wu et al., 2016), which may themselves be biased.In particular, in Utah and Nevada, simulated near-surface dust concentrations are about 2-3 times as large as observed, indicating significant overestimation of dust emissions in the region.

Aerosol-in-snow concentrations
Figure 4 shows the spatial distributions of BC and dust mass mixing ratios in snow in winter and spring from VR-CESM simulations.The BC-in-snow mass mixing ratio in the Rocky Mountains ranges from 2-50 ng g −1 , which is consistent with a previous study (Qian et al., 2009).The dust-insnow mass mixing ratio (0.1-50 µg g −1 ) is about 2-3 orders of magnitude higher than that of BC-in-snow.The spatial pattern of BC-in-snow mixing ratios is consistent with that of the near-surface atmospheric BC concentration, which features higher values in northern Utah and southern Idaho and lower values in the higher mountains (Fig. 2b).The dust-insnow mixing ratios are higher in Utah and downwind regions (western Colorado and southern Idaho), which is consistent with the distribution of near-surface atmospheric dust con- centrations.Dust-in-snow mixing ratio is also higher in the northern Great Plains, where dust emission is also evident (Fig. 2c).In addition, BC and dust mixing ratios are larger (10-100 ng g −1 and 2-50 µg g −1 , respectively) in the Southern Rockies than in the Northern Rockies and Greater Yellowstone region.BC and dust mass mixing ratios are smaller in the Greater Yellowstone region with ranges from 10-50 ng g −1 and from 0.2-2 µg g −1 , respectively, and are smallest in the Northern Rockies with values below 20 ng g −1 and below 2 µg g −1 , respectively.BC and dust mixing ratios in snow are larger in spring than in winter in most of the Rocky Mountain region.This is due to larger deposition of BC and dust in spring than in winter, resulting from larger northward transport of BC and dust in spring (Fig. not shown).Larger dust deposition in spring can also be partly explained by the larger dust emission in this season.
The comparison of BC mass mixing ratios in the snow column at the 17 sites from VR-CESM simulations and observations is shown in Fig. 5. Observed BC-in-snow mass mixing ratios range from 5.5 to 33.6 ng g −1 at the 17 sites.Simulated BC mixing ratios range from 8.3 to 30.6 ng g −1 at these sites, which are in the range of observations.Despite this, simulated BC-in-snow mass mixing ratios differ from the observations by a factor of up to 4 at some stations.Averaged across all the 17 sites, the simulated BC mass mixing ratio is 35 % larger than the observed value.Note that there is a large interannual variability of BC-in-snow mass mixing ratios, such as at site 16, as shown in Doherty et al. (2016).Therefore, the observation period was not long enough for the derivation of a climatological mean as in the simulation, which may partly explain the inconsistency between the observations and simulations.
Note that although near-surface BC concentrations in the atmosphere are underestimated in the Rocky Mountains in the model (Sect.4.1), BC mass mixing ratios in the snow are not overall underestimated.In our study, we calibrate the observational data of Doherty et al. (2014) by using a correction factor based on the comparison of ISSW and SP2, assuming that SP2 can more accurately measure the mass mixing ratio of BC compared to ISSW (Doherty et al., 2016).However, although SP2 can provide a direct measurement of BC, SP2 may underestimate the real amount of BC-in-snow mass when BC is attached to larger particles (e.g., dust and sea salt) or aggregates to large sizes in snow due to the size range (e.g., ∼ 0.08-0.7 µm) limitations in SP2 (Qian et al., 2015).Because of this, the real amount of BC-in-snow mass may be higher than that measured by SP2.Another reason for the inconsistency of BC mass mixing ratios in snow and near-surface BC concentrations in the atmosphere may be related to the compensating errors in BC deposition and snowfall.This inconsistency may also be related to snow aging and melting and BC-in-snow accumulation and flushing out, which are associated with large uncertainties (Flanner et al., 2007;Qian et al., 2014).
For dust-in-snow, the simulated mean dust mass mixing ratio in snow in May-June is 31.0(27.8) µg g −1 in the San Juan Mountains (Grand Mesa), with the standard deviation, minimum, and maximum being 20.4 (10.0), 8.9 (14.8), and 81.4 (50.4) µg g −1 , respectively.These values are 1 to 2 orders of magnitude smaller compared to the observed mixing ratios from Skiles et al. (2015), which showed that, at the end of the snow season, the total dust-in-snow mass mixing ratios range from 0.2 to 4.8 and from 0.6 to 1.7 mg g −1 , respectively, in the San Juan Mountains and the Grand Mesa.Much smaller dust-in-snow mass mixing ratios in the simulations may be ascribed to the fact that the model only accounts for dust particles with diameters smaller than 10 µm, while the observations include all sizes of dust particles in the snow.Observation by Reynolds et al. (2016) in the Colorado region showed that the mass concentrations of dust particles in snow are mostly from larger particles with diameters larger than 10 µm.Therefore, the model may underestimate the impacts of dust deposition into snow.The dust impacts calculated in this study, which will be discussed below, should be regarded as those from dust particles with diameters smaller than 10 µm.

Surface radiative effect (SRE) by aerosol-in-snow
Figure 6 shows the spatial distribution of the instantaneous surface radiative effect (SRE) due to BC-and dust-induced snow albedo change, respectively, in winter (December-January-February) and spring (March-April-May).Due to the decrease in surface albedo, surface net shortwave radiation is increased.The spatial patterns of SRE are determined by both the amount of aerosol-in-snow and the snowpack distribution (snow depth and snow cover fraction).Finerscale structures of SRE in the Rocky Mountains and the adjacent regions are simulated by VR-CESM with a higher horizontal resolution compared to previous simulations by coarse-resolution GCMs (e.g., Flanner et al., 2009;Yasunari et al., 2015).The SRE is generally above 0.2 W m −2 over the mountains, especially in the Greater Yellowstone region and Southern Rockies.SRE can reach similar magnitudes on the southern periphery of the Northern Rockies and the western side of the Greater Yellowstone region, where higher near-surface atmospheric concentrations and in-snow mass mixing ratios of BC and dust are simulated (Figs. 2  and 4).SRE is stronger in spring than in winter for both BC and dust, which is consistent with previous studies (Flanner et al., 2009;Yasunari et al., 2015).This is because of the stronger solar insolation and larger albedo reduction due to snow aging, aerosol accumulation within snow, and feed- backs in spring.Dust emissions, and consequent dust transport and deposition, are higher in spring than in winter, which may also partly contribute to the larger dust-induced SRE in spring than in winter.BC-induced SRE is somewhat larger than dust-induced SRE in both winter and spring.BCinduced SRE is mostly below 1 W m −2 in winter, but reaches up to 2-5 W m −2 in spring.Dust-induced SRE is mostly below 0.5 W m −2 in winter and increases to 1-5 W m −2 in spring.
Compared to the Greater Yellowstone region and Southern Rockies, SRE in the Northern Rockies is much smaller (mostly below 0.05 and 0.5 W m −2 in winter and spring) because of smaller aerosol-in-snow mixing ratios in this region (Fig. 4).Note that BC-induced SRE is still significant (mostly around 0.2-2 and 2-5 W m −2 in some local regions) in the northern Great Plains, eastern US, southern Canada, and eastern Canada.This was also shown in previous studies using coarse-resolution GCMs (e.g., Flanner et al., 2009;Yasunari et al., 2015).In addition, the model also simulates non-negligible dust-induced SRE (mostly around 0.05-0.2and up to 0.2-0.5 W m −2 in some local regions) near the dust sources in southern Canada and the northern Great Plains.
Figure 7 shows the monthly variations in SRE induced by BC and dust SDE in the five regions (Northern Rockies, Greater Yellowstone region, Southern Rockies, eastern Snake River Plain, and southwestern Wyoming).Table 2 gives the regionally averaged winter and spring SRE in these five regions.Consistent with the spatial distributions shown in Fig. 6, aerosol-induced SRE averaged in the Northern Rockies is about one-half to one-fourth of that in the Greater Yellowstone region and Southern Rockies.Compared to that in winter, SRE is much larger in spring, which is a result of aerosol accumulation in snow and relatively strong solar insolation.Maxima in the monthly SRE occur in April-May in the three mountainous regions (Northern Rockies, Greater Yellowstone region, and Southern Rockies), which is consistent with the progress of snowmelt after the peaks of snow water equivalent (early to middle April; see Fig. 11 of Wu et al., 2017).In the eastern Snake River Plain and southwestern Wyoming, maxima in monthly SRE occur in March, which is different from the three mountainous regions because the snowmelt period begins earlier (in February to March) in these two regions (Sect.4.4).Regional mean total SRE in spring induced by BC and dust can reach up to 1.58-1.7 W m 2 with monthly peaks around 2.0 W m −2 in the Greater Yellowstone region and Southern Rockies.In the eastern Snake River Plain and southwestern Wyoming regions, regional mean total SRE in winter and spring is around 0.4-0.6 and 0.7-0.8W m −2 , respectively.Dust-induced springtime SRE can contribute to about 20-30 % of total springtime SRE in the northern part of the Rocky Mountains (Northern Rockies, Greater Yellowstone region, and eastern Snake River Plain).In the southern part of Rocky Mountains (southwestern Wyoming and Southern Rockies), dust-induced springtime SRE contributes more significantly (about 30-40 %) to total springtime SRE.Note that the dust-induced SRE shown here does not take into account dust particles larger than 10 µm, which may constitute the majority of dust-in-snow mass (Reynolds et al., 2016).Therefore, our estimations of dust-induced SRE may be biased low.

Impacts of aerosol SDE on the surface temperature and snowpack
Figure 8 shows surface air temperature, snow water equivalent, and snow cover fraction changes due to the aerosol SDE in winter and spring.Snow water equivalent is defined as the amount of water contained within the snowpack measured in kg m −2 , which is equivalent to millimeters after dividing by the density of water (1000 kg m −3 ).The snow cover fraction is defined as the fraction of surface area covered by snow.These changes are derived from the difference between the two simulations (CTL and NoSDE).The crosses in the figure denote the regions where changes are statistically significant at the 0.1 level.Although SRE is largest over the mountains, surface air temperature change is largest in the regions adjacent to the mountains, such as over the eastern Snake River Plain, northern Utah, and central and southwestern Wyoming, where surface air temperatures are increased by around 0.5-2 • C due to the aerosol SDE.The large surface air temperature increase corresponds well with the significant reductions of snow water equivalent (by 2-50 mm) and snow cover fraction (by 5-20 %) in these regions.This indicates a pronounced positive feedback between snow albedo, radiation, and surface temperature in the regions adjacent to the mountains, where snow water equivalent values are relatively lower and snow cover fractions are smaller than those over the mountains.The positive feedback amplifies the surface warming and snow melting, as was also found in a previous study using the Weather Research and Forecasting (WRF) model (Qian et al., 2009).We note, however, that both snow water equivalent and snow cover fraction are larger over the mountains.For example, winter and spring snow water equivalent is mostly above 50 mm on the high mountains (see Fig. 8 of Wu et al., 2017).Local aerosol SDE may also induce substantial impacts on the surface temperature and snowpack on the high mountains, but these impacts may be canceled out by the increase in snowfall (Fig. 9).As shown in Fig. 8, the smaller change in surface air temperature over the mountains corresponds well with the increase in snow water equivalent and snow cover fraction (especially in the Northern Rockies and the Greater Yellowstone region).
The increase in snowfall in Fig. 9 is likely related to the large-scale circulation change due to aerosol SDE. Figure 10 shows wintertime tropospheric temperature and zonal winds in CTL and NoSDE simulations and their difference.In the NoSDE simulation, we have turned off the SDE not only in the Rocky Mountain region, but also in other regions of the globe.Due to aerosol SDE, temperature is increased in the high latitudes of the Northern Hemisphere (Fig. 10c), which can reduce the meridional temperature gradient, thus leading to a weakening polar jet stream north of 50 • N (Fig. 10f).This suggests a shift to a more meridional wind pattern in winter, which can enhance the broader meanders and thus the formation of winter storms (Wu et al., 2017).Enhanced winter storm activity further reduces surface temperature to the north of the Rocky Mountain region and in the northern part of the Rocky Mountain region (Figs.8a and 10c).This, together with the increased temperature in the southwestern US and southern part of the Rocky Mountain region, increases the meridional temperature gradient and leads to a stronger westerly at 30-45 • N (Fig. 10f).A stronger westerly at 30-45 • N favors water vapor transport from the Pacific Ocean.The enhancement of winter storm activity and water vapor transport may lead to an increase in precipitation (mainly in terms of snowfall in winter).In spring, the change in temperature and zonal winds is similar to that in winter, but with a northward shift of the patterns as a result of northward movements of the polar jet stream and westerlies in spring (not shown).Therefore, the change in snowfall is likely a result of circulation change induced by SDE from both the Rocky Mountain region and remote regions.It is worth isolating the impacts of SDE from the Rocky Mountain region and remote regions (e.g., high latitudes) in the future.Note that increases in snow water equivalent and snow cover fraction in the Northern Rockies and Greater Yellowstone region due to aerosol SDE do not pass the significance test at the 0.1 level because of the large interannual variability in these regions.
Table 2 gives the winter and spring surface air temperature changes due to LAAs in snow averaged over the five regions.The seasonal mean surface air temperature change is around 0.9-1.1 • C in the eastern Snake River Plain in winter and spring, while this change is around 0-0.2 • C (winter) and around 0.3-0.5 • C (spring) in the mountainous regions (Northern Rockies, Greater Yellowstone region, and Southern Rockies).In Table 2, we also show the efficacy of snow albedo forcing, which is defined as the ratio of local surface air temperature change to SRE over a specific region.The efficacy is mostly around 0.1-0.6 in the three mountainous regions, but it is 1.3-2.2 in the eastern Snake River Plain and southwestern Wyoming.This indicates that stronger snow albedo feedbacks exist in the latter two regions.
Figures 11-12 show the monthly evolution of regional mean surface air temperature, snow water equivalent, snow cover fraction, and their changes due to aerosol SDE in the eastern Snake River Plain and southwestern Wyoming.Monthly variations in surface air temperature, snow water equivalent, and snow cover fraction are similar between the two regions: lowest surface air temperature and largest snow cover fraction in January and highest snow water equivalent in February-March.Significant changes in these variables from the aerosol SDE occur in both regions.The largest surface air temperature increase is 1.5 • C in the eastern Snake River Plain and 1.6 • C in southwestern Wyoming, occurring in April and December, respectively.In the eastern Snake River Plain (southwestern Wyoming), aerosol SDE leads to the reduction of snow water equivalent by 6-28 mm (6-16 mm) and snow cover fraction by 5-15 % (6-19 %) from December to March.In April (late snowmelt period) when snow water equivalent and snow cover fraction are both relatively small, the aerosol SDE is more significant, which reduces snow water equivalent values and snow cover fractions by about half.

Runoff change induced by aerosol SDE
In the model, the total runoff includes the surface runoff and subsurface runoff.Our simulations show that the spatial distribution and seasonal evolution of surface runoff and subsurface runoff is generally similar to total runoff (not shown), and surface runoff and subsurface runoff account for 30- 40 and 60-70 %, respectively, of annual total runoff in the mountains.Here we only show the simulation results of total runoff, as both surface runoff and subsurface runoff will flow to rivers and become discharge.Runoff is mainly from rainfall and snowmelt.The change in rainfall is shown in Fig. 9cd, and the snowmelt change is shown in Fig. 13.Aerosol SDE increases the snowmelt by 0.1-2 mm day −1 in the mountains during the snow accumulation and early snowmelt period (in autumn, winter, and spring).In the late snowmelt period, aerosol SDE reduces the snowmelt due to less snowpack available for melting in the plains (in spring) and in the mountains (in summer).Note that snowmelt is slightly reduced by the aerosol SDE in autumn in the Southern Rockies, which is a result of less snowpack available for melting due to the reduced snowfall in this region (not shown).
Because of the change in rainfall and snowmelt due to aerosol SDE, runoff changes too. Figure 14 shows the runoff change induced by aerosol SDE in four seasons.In winter, runoff is barely modified by the aerosol SDE in the Rocky Mountains, except in the Northern Rockies where runoff is increased by 0.1-2 mm day −1 and associated with increased rainfall (Fig. 9c) and increased snowmelt (Fig. 13a).In spring, runoff changes the most compared to all of the other seasons, with the runoff increased by up to 0.5-2 mm day −1 in the mountainous regions.This is mainly due to the increase in snowmelt resulting from surface warming (Fig. 13b) and due to more snow available for melt resulting from snowfall increase (Fig. 9f).The changes in runoff are statistically significant at the 0.1 level in most of the mountainous regions in spring.Absolute runoff increases are stronger in the Northern Rockies and Greater Yellowstone regions than in the Southern Rockies in terms of the area and magnitude, probably due to the smaller snow water equivalent in the Southern Rockies (Wu et al., 2017).As more snowmelt occurs in spring, less snowpack is available for melt in summer, and thus surface runoff is reduced by about 0.1-1 mm day −1 .There is little runoff change in autumn, as there is less runoff generated from rainfall and snowmelt than in other seasons.Overall, BC and dust residing in snow accelerate the hydrologic cycles by increasing the runoff in spring and reducing the runoff in summer.Surface warming also increases the ratio of rainfall to total precipitation, which can accelerate the generation of runoff.Note that in some regions of the plains, such as central and eastern Montana, southwestern Wyoming, and the Snake River Plain, the snowmelt changes by 0.1-1 mm day −1 due to the aerosol SDE, but the runoff changes little.This is because the water generated from snowmelt is mainly stored in soil or transformed into evapotranspiration in these regions.Also note that there are statistically significant increases in runoff in the southern Great Plains in spring, but the change is small (around 0.05-0.1 mm day −1 ; Fig. 14b).This change is a result of slight increases in both rainfall (Fig. 9d) and snowmelt (Fig. 13b).
Figure 15 shows the monthly evolution of runoff and its change due to the aerosol SDE in the three mountainous re- gions (Northern Rockies, Greater Yellowstone region, and Southern Rockies).In the three regions, runoff peaks in the late spring and early summer (in May in the Northern Rockies and Southern Rockies and in June in the Greater Yellowstone region) when snow melting progresses after the peak of snow water equivalent in early to middle April (Wu et al., 2017).This indicates the significant contribution of snowmelt to runoff.Overall, runoff changes are larger in the Northern Rockies and Greater Yellowstone region than in the Southern Rockies, which is consistent with the spatial distribution of runoff changes shown in Fig. 14.Runoff is significantly increased in spring and decreased in June and July, indicating the acceleration of the local hydrologic cycle by aerosol SDE.In the Northern Rockies, runoff is also increased from October to March but in much smaller magnitudes (below 0.2 mm day −1 ) compared to April and May.In April (May), runoff is increased by 0.39 (0.56), 0.22 (1.00), and 0.17 (0.15) mm day −1 in the Northern Rockies, Greater Yellowstone region, and Southern Rockies, respectively.This increase contributes to 26 (13 %), 42 (27 %), and 29 % (7 %) of the runoff from the NoSDE simulation in April (May) for the three regions, respectively.The reduction of runoff in June is relatively small (0.06 and 0.11 mm day −1 , respectively) in the Northern Rockies and Greater Yellowstone, only accounting for 2 % of runoff from the NoSDE simulation.However, it reaches up to 0.18 mm day −1 in the Southern Rockies, which accounts for 15 % of runoff.In addition, due to the reduction of snow available for melting later in July, the runoff is further reduced.Runoff is relatively smaller in July versus previous months, and aerosol SDE can reduce the runoff by 0.04 (8 %), 0.17 (23 %), and 0.06 mm day −1 (16 %) in the three regions, respectively.Note that due to the increase in precipitation, the annual mean runoff is increased by 0.12 (12 %), 0.09 (10 %), and 0.01 mm day −1 (2 %) in these three regions, respectively.

Conclusions
In this study, we use VR-CESM to quantify the impacts of LAA (BC and dust) deposition to the snowpack, hydrologic cycles, and surface air temperatures in the Rocky Mountains.Our previous study has shown that VR-CESM reasonably reproduces the spatial distributions and seasonal evolution of snowpack in the Rocky Mountains (Wu et al., 2017).Here we show that the model simulates similar magnitudes of near- Due to the deposition of LAAs to snow, surface net shortwave radiation is increased.Regionally and seasonally averaged SRE induced by LAAs in snow is 0.1-0.5 W m −2 in winter in the three mountainous regions (Northern Rockies, Greater Yellowstone region, and Southern Rockies) and 0.4-0.6W m −2 in the two regions around the mountains (eastern Snake River Plain and southwestern Wyoming).Seasonal average SRE is much larger in spring and reaches up to 0.6-1.7 W m −2 in these five regions (Table 2).Dust contributes 21-42 % to the total SRE induced by LAAs in snow in spring, indicating the important role of dust residing in snow.Of the five regions, dust contributes the most (42 %) to the total SRE in the Southern Rockies.This is not unexpected as this region is close to dust sources in the Colorado Plateau.
As a result of SRE induced by LAAs in snow, surface air temperature increases in most of the Rocky Mountain region.The surface air temperature increase is largest over the eastern Snake River Plain and southwestern Wyoming, with winter and spring surface air temperature increased by 0.9-1.1 • C. Significant reductions of snow water equivalent (by 2-50 mm) and snow cover fraction (by 5-20 %) occur in these two regions, indicating a strong positive snow albedo feedback there.
Aerosol SDE accelerates the hydrologic cycle in the mountainous regions.In April and May, monthly mean runoff is increased by 7-42 % in the three mountainous regions (Northern Rockies, Greater Yellowstone region, Southern Rockies).This is because of the accelerated snowmelt resulting from surface warming and the increased snowfall resulting from enhanced winter storm activity and water vapor transport from the Pacific Ocean.This enhancement may be related to large-scale circulation changes.In the later stage of snowmelt, monthly runoff is reduced by 2-15 % in June and 8-23 % in July in the three mountainous regions.In particular, aerosol SDE leads to a reduction of total runoff by about 15 % in June and July in the Southern Rockies.This highlights the important role of aerosol SDE in modulating the hydrologic cycle in these mountainous regions.
We note that VR-CESM still underestimates the nearsurface BC concentrations; however, it overestimates BC-insnow concentrations by 35 % for the average across the 17 observational sites.For dust-in-snow, the model used in this study only accounts for dust particles smaller than 10 µm, while observations made by Reynolds et al. (2016) suggest that most airborne and in-snow dust mass concentrations are characterized by dust particles with diameters larger than 10 µm in the Utah-Colorado region.Therefore, our simulations may significantly underestimate the impacts of dust-insnow, especially over the Southern Rockies.In the Southern Rockies, our simulations suggest that SRE induced by dust-in-snow can reach up to 2-5 W m −2 , which is nearly an order of magnitude smaller than values given by Painter et al. (2007) and Skiles et al. (2015) based on observed dustin-snow particles in the same region.Note that such bias in SRE may become smaller in the Greater Yellowstone region and Northern Rockies as these regions are farther from the dust source regions than the Southern Rockies.Future observations of LAAs in snow, particularly for the temporal evolution of LAAs in different snow layers, and detailed size distribution measurements of dust particles in snow will help reduce the uncertainties in the model quantification of the impacts of LAAs in snow.
Although uncertainties still exist, our results show that LAAs in snow can significantly affect the snowpack and consequent hydrologic cycle in the Rocky Mountains.Previous studies have demonstrated that snowpack on the Rocky Mountains has declined significantly in the second half of 20th century (e.g., Pederson et al., 2011).The role of LAAs in this decrease in snowpack is still unknown.It would be interesting to investigate the role of LAAs and compare it with those of other climate factors (such as natural climate variability and greenhouse gas concentrations).Moreover, BC and dust emissions may also be subject to changes in the future.Therefore, for better projections of future changes in Rocky Mountain snowpack, the impacts of LAAs in snow under future emission scenarios need to be taken into account.
of the Meteorological Public Welfare Profession of China (grant GYHY01406021), the National Natural Science Foundation of China (grant 41575095), and the Chinese Academy of Sciences "The Belt and Road Initiatives" program on international cooperation: Climate Change Research and Observation Project (grant 134111KYSB20160010).We thank the team for maintaining the Interagency Monitoring of PROtected Visual Environments (IMPROVE) network and making the observation dataset available to use (http://vista.cira.colostate.edu/Improve/improve-data/).We also thank Sarah Doherty from the University of Washington and Sara McKenzie Skiles from the University of Utah for providing the observations of absorbing aerosols in snow and helpful suggestions on the use of the data.We thank Alan M. Rhoades and Paul A. Ullrich from the University of California, Davis and Colin M. Zarzycki from NCAR for helpful discussions during this study.We would like to acknowledge the use of computational resources for conducting the model simulations (ark:/85065/d7wd3xhc) at the NCAR Wyoming Supercomputing Center provided by the NSF and the state of Wyoming and supported by NCAR's Computational and Information Systems Laboratory.
Edited by: Yafang Cheng Reviewed by: two anonymous referees

Figure 1 .
Figure 1.(a) Model meshes for variable resolution (uniform 1 • with refined 0.125 • in the Rocky Mountains) used in VR-CESM.Note that each element shown contains additional 3 × 3 collocation grid cells.(b) Terrain height (m) in the western US for the variable-resolution grid used in VR-CESM (Sect.2).The refined region at a resolution of 0.125 • is surrounded by dashed lines.(c) Five regions are identified for the analysis in this study, including three mountainous region (1, Northern Rockies; 2, Greater Yellowstone region; 3, Southern Rockies) and two regions in the plains around the mountains (4, eastern Snake River Plain; 5, southwestern Wyoming).

Figure 2 .
Figure 2. Spatial distribution of cold season (winter and spring) mean (a) BC emission flux and (b) near-surface BC concentration from the VR-CESM simulation; (c) and (d) for dust emission flux and near-surface dust concentration, respectively.Also shown are the IMPROVE stations (blue open circle) selected for model validation, with the size of the circles from small to large indicating the magnitude of observed near-surface concentrations of BC (b) and dust (d).The black rectangles in (b) and (d) denote the five regions (A, West Coast; B, Rocky Mountains; C, Utah and Nevada; D, Southwestern US; E, Great Plains), which will be used to classify the stations in Fig. 3.Note that units for BC and dust concentrations are ng m −3 and µg m −3 , respectively.

CFigure 3 .
Figure3.Comparison of cold season (winter and spring) mean near-surface (a) BC and (b) dust concentrations at IMPROVE stations from VR-CESM simulation and IMPROVE observations.Also given are the mean results at all the stations from simulation and observations and their correlation coefficient (R).The 1 : 1 (solid) and 1 : 5/5 : 1 (dashed) lines are plotted for reference.

°°°°°°°°°°°°°°°°°°°°F igure 4 .
Winter (December-January-February, DJF; a, c) and spring (March-April-May, MAM; b, d) mean BC (a, b) and dust (c, d) mass mixing ratios in the snow column.Also shown are the stations for observations of BC mass in the snow column (a, b) and for observations of dust mass in the snow column in the San Juan Mountains and the Grand Mesa (denoted by the diamond and square, respectively; c, d).

Figure 5 .
Figure 5.Comparison of BC mass concentrations in the snow column (C BC ) at the 17 sites (seeTable 1) from VR-CESM simulations and observations, with the error bars denoting the corresponding standard deviations.The observations are compiled from previously published studies (Table1).If multiple observations are recorded at a certain site, the observed standard deviations are calculated from these multiple observations (Sect.3).Simulated BC mass concentration in the snow column and its standard deviation are calculated from the 25-year mean and standard deviation of simulation on the same month and day as the observations (Sect.3).The 1 : 1 (solid) and 1 : 5/5 : 1 (dash) lines are plotted for reference.

Figure 7 .
Figure7.Monthly variations in surface radiative effect (SRE; W m −2 ) during the water year (1 October to 30 September) averaged over the Northern Rockies, Greater Yellowstone region, Southern Rockies, eastern Snake River Plain, and southwestern Wyoming.

Figure 8 .
Figure 8. Changes in surface air temperature (a, b; • C), snow water equivalent (c, d; mm), and snow cover fraction (e, f; %) in winter (a, c, e) and spring (b, d, f) induced by BC-in-snow and dust-in-snow.The crosses denote the regions where changes are statistically significant at the 0.1 level.

Figure 10 .
Figure 10.(a-c) Wintertime temperature ( • C) and (e-f) zonal winds (m s −1 ) averaged at 102-125 • W from CTL and NoSDE simulations and their difference.Note that zonal winds are averaged for a range of longitudes, which correspond to the east-west boundary of the western US, including the Rocky Mountains and upwind regions (Fig. 1).

Figure 11 .
Figure 11.Seasonal evolution of (a) surface air temperature, (c) snow water equivalent, and (e) snow cover fraction and their changes due to SDE (b, d, f) averaged over the eastern Snake River Plain.

Figure 13 .
Figure 13.Snowmelt change (mm day −1 ) due to SDE of BC and dust in four seasons: (a) December-January-February (DJF), (b) March-April-May (MAM), (c) June-July-August (JJA), and (d) September-October-November (SON).The crosses denote the regions where changes induced by SDE are statistically significant at the 0.1 level.

Figure 15 .
Figure 15.Seasonal evolution of total runoff including surface and subsurface runoff (a, c, e) and their change (b, d, e) in the Northern Rockies (a, b), the Greater Yellowstone region (c, d), and the Southern Rockies (e, f).The unit is mm day −1 .

Table 1 .
Observations of BC mass concentration in the snow column (C BC , ng g −1 , i.e., ng gram BC per g snow) in the Rocky Mountain region compiled from previously published literature.

Table 2 .
Skiles and Painter (2016b)y et al. ( , 2016) )(March-April-May) mean surface shortwave radiative effect (SRE; W m −2 ) due to BC alone, dust alone, and BC and dust together in snow, as well as surface air temperature (SAT; • C) change and the efficacy of SRE in SAT change in the five regions (see Fig.1c).Note that SRE induced by BC and dust together is slightly larger than the sum of SRE induced by BC and SRE by dust separately.BC by ISSW to C BC by SP2 based on their linear relationship for the estimation of real C BC .This calibration is applied to the dataset ofDoherty et al. (2014)in our study, and thus the observations ofDoherty et al. (2014Doherty et al. ( , 2016) )used here are comparable.In addition,Skiles and Painter (2016b)made daily measurements of BC-in-snow with an SP2 in the Senator Beck Basin Study Area (SBBSA) in the San Juan Mountains during a period of 2 months (late March to middle May) in 2013.
Doherty et al. (2016) by BC (dust) to the sum of SRE by BC and SRE by dust is given in parentheses next to SRE by BC (dust).bTheefficacy of snow albedo forcing ( • C increase per 1 W m −2 ) is defined as the ratio of SAT change to SRE. variations in C BC at four stations, three in Idaho (January to March of 2014) and one in Utah (January to February of 2013 and 2014).Doherty et al. (2016)also calibrated the ISSW measurements using an incandescence technique (a singleparticle soot photometer, SP2) in a subset of the observations, which was intended to capture C BC more accurately, and derived a ratio of C surface dust concentrations at most stations in the Rocky Mountain region compared to IMPROVE observations.The model tends to underestimate near-surface atmospheric BC concentrations mostly by a factor of 1.5-5 in the Rocky Mountain region.The underestimation of near-surface BC concentrations may be due to the absence of local sources in the BC emissions dataset used and too-weak transport in the model.Simulated aerosol-in-snow concentrations are closely related to the distributions of both snowpack and near-surface atmospheric aerosol concentrations.Simulated BC-in-snow concentrations range from 2 to 50 ng g −1 in the Rocky Mountain region, and they are 35 % larger than the observations for the average at the 17 sites.