Examining the atmospheric radiative and snow-darkening effects of black carbon and dust across the Rocky Mountains of the United States using WRF-Chem

Affiliations 1Institue of the Environment and Sustainability, University of California Los Angeles, Los Angeles, California, 90095, U.S.A. 2Department of Atmospheric Science, University of Wyoming, Laramie, Wyoming, 82071, U.S.A. 10 3Department of Atmospheric Sciences, Texas A&M University, College Station, Texas, 77843, U.S.A. 4School of Earth and Space Sciences, University of Science and Technology of China, Hefei, China 5Anhui Province Key Laboratory of Extreme Events in Atmosphere, University of Science and Technology of China, Hefei, China 15 Correspondence to: Xiaohong Liu (xiaohong.liu@tamu.edu)


Introduction
The arid Rocky Mountains of the western U.S. (WUS) receive most of their precipitation in the form of snowfall from October through March. The resulting snowpack forms a network of natural mesoscale storage reservoirs that provide approximately 80% of the surface water across the region during 35 the warm season (Serreze et al., 1999;Hamlet et al., 2007). All life in the region fundamentally depends on the timed release of runoff from snowpack; humans rely on these resources for agriculture and power generation. In recent decades however, there have been observed changes in the hydrology across the WUS associated with external climate forcings (e.g., anthropogenic climate change) that may be acting to compromise the security of water resources across the region and beyond (Hamlet et al., 2007). 40 Numerous studies have shown that annual maximum snow water equivalent (SWE) values have decreased since 1950 (Das et al., 2009;Mote 2006;Pierce et al., 2008), increasing (decreasing) runoff discharge in the winter and spring (summer) (Rajagopalan et al., 2009;Qian et al., 2009). Externally forced warming associated with greenhouse gases and light-absorbing aerosols (LAA; Pierce et al., 2008) and LAA deposition on snowpack (Flanner et al., 2007;Qian et al., 2009, rather than natural 45 climate variability, are believed to be the major contributors to this decrease. LAA such as black carbon (BC) and dust (refereed to collectively as BCD) can affect the hydrology across the WUS as they interact with sunlight, altering the vertical thermodynamic structure of the atmosphere. These aerosol-radiation interactions (ARI) may lead to changes in precipitation amount and type (Pederson et al., 2011). BCD may also deposit on snowpack, increasing the surface absorptivity 50 and enhancing melting in a process referred to as the snow darkening effect (SDE; Warren and Wiscombe, 1980;Painter et al., 2007). Surface warming is generally brought about by the SDE, while the surface can either cool or warm from ARI. Both effects have been shown to be important across the region however, especially since the surface radiative budget is sensitive to small perturbations in albedo (Pepin et al., 2015). 55 BCD aerosols find their way into the WUS from both near-field and far-field sources. BC, produced via the incomplete combustion of fossil fuels and biomass burning, is primarily emitted in WUS cities and by wildfires (Bond et al., 2013). Dust meanwhile is primarily emitted from southwestern U.S. deserts via wind erosion (Tegen et al., 2004). Following emission, both aerosols are transported https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. downstream by prevailing westerlies toward the Rocky Mountain region. Additionally, a sizable 60 component of atmospheric dust across the WUS originates from Asian sources (Fischer et al., 2009). BCD SDEs and ARI across the WUS have been studied using global climate models (GCMs) (Flanner et al., 2007;Qian et al., 2009; and regional climate models (RCMs) (Wu L. et al., 2018), each with their own advantages and drawbacks. Heterogeneous mesoscale meteorology features (i.e. precipitation, temperature, and snow characteristics) can be simulated better with RCMs than GCMs, 65 as higher grid resolutions are typically used. Wet removal by precipitation rather than dry removal is a more effective pathway for cleansing the atmosphere of BCD (Zhao et al., 2014). Therefore, highresolution (cloud-resolving) simulations, through their more physically based treatment of orographically forced precipitation, should better simulate these aerosols' lifecycle. Additionally, RCMs such as the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem; Grell et al. (2005)) have 70 chemical and aerosol options that are generally more sophisticated than typical GCMs.
While WRF-Chem has been used to study the impacts of BCD SDEs and ARI across California on convective-permitting scales (Wu L. et al., 2018), the application of this model to the American Rockies as a whole has not been made. Smaller inner-continental cities and municipalities across this zone may be highly sensitive to changes in hydrology brought about by BCD effects, hence providing motivation for this 75 study.
Using WRF-Chem, we seek to quantify perturbations to WUS meteorology and hydrology induced by BCD for water year 2009. Sensitivity experiments are run on convective-allowing scales (4 km grid spacing) to isolate the effects of BCD SDEs and ARI on temperature, precipitation, snow, and runoff.
This study begins with a presentation of the model, methodology, and data in section 2. Section 3 80 provides a meteorological and chemical validation of WRF-Chem. The radiative effects of BCD associated with SDEs and ARI are explored in section 4, and their effects on WUS weather are examined in section 5.
Section 6 briefly evaluates the implications of undersimulated dust emissions, while Section 7 examines the combined effects (SDE+ARI) of BCD. Conclusions are presented in Section 8. In this study, WRF-Chem 3.5.1 updated by University of Science and Technology of China (USTC) is used. This USTC version of WRF-Chem includes capabilities not available in publicly released WRF-Chem versions such as the diagnosis of radiative effects of aerosol species, land surface coupled 90 biogenic VOC emissions, and aerosol-snow interactions (Zhao et al., 2013a(Zhao et al., ,b, 2014(Zhao et al., , 2016Hu et al., 2019).) The model is run on a 4-km grid with 50 vertical levels across a large portion of the WUS ( Figure   1). The SNow, ICe, and Aerosol Radiative (SNICAR) model (Flanner et al., 2007), which uses BCD deposition flux from the atmosphere to treat the SDE and treats ice and snowpack aging, was coupled into the Community Land Model version 4.0 (CLM4; Oleson et al., 2010 ) by Zhao et al. (2014). SNICAR uses 95 snow water (both ice and liquid) and aerosol loading information to compute the snowpack's radiative properties within multiple snow layers based on the theories by Warren and Wiscombe (1980) and Toon et al. (1989). The utility of SNICAR in simulating the albedo reductions in snow has been tested in laboratory experiments (Hadley and Kirchstetter, 2012).
The MOdel for Simulating Aerosol Interactions and Chemistry (MOSAIC; Zaveri et al., 2008) and 100 the Carbon Bond Mechanism version "Z" (CBM-Z; Zaveri and Peters, 1999) photochemical model are used to treat aerosol and atmospheric chemistry. MOSAIC uses a 4-bin sectional approach to simulate the aerosol size distributions of BC, dust, sulfate, ammonium, nitrate, organic matter, and sea salt for radii of 0.039-10 m. Additionally, MOSAIC treats the processes of aerosol nucleation, coagulation, condensation, water uptake, and aqueous chemistry. Aerosol dry deposition is handled via the method in Binkowski and 105 Shankar (1995), which includes Brownian and turbulent diffusion as well as gravitational settling. Wet deposition of aerosols and gases by in-cloud and below-cloud precipitation scavenging is treated following Easter et al. (2004). Similar to Zhao et al. (2013a) and Wu et al. (2017;, aerosols are assumed to be internally mixed within each size bin. Aerosol optical properties such as extinction, single-scatter albedo, and asymmetry parameter are computed as a function of wavelength at each grid point based on the bin-110 and volume-averaged refractive index for each aerosol species (Fast et al. 2006). ARI are treated in the radiation code via the methodologies in Zhao et al. (2011), in which the direct radiative effects is computed diagnostically using the method in Ghan et al. (2012). Aerosol radiative feedbacks are enabled, and the number concentration of activated aerosols is treated prognostically in the cloud microphysics scheme (Morrison et al. 2009) following Gustafson et al. (2006). 115 https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License.

Emission data
Anthropogenic emissions from the Environmental Protection Agency's (EPA) 2011 National Emission's Inventory (EPA NEI-11; https://www.epa.gov/air-emissions-inventories/2011-nationalemissions-inventory-nei-data) are used. These emissions contain location-specific point and area source 120 emissions and are interpolated to a 4-km grid using the open-source software emiss_v04.F (ftp://aftp.fsl.noaa.gov). Biomass burning emissions, available on a ~1 km grid from the Fire INventory from the National Center for Atmospheric Research (NCAR) (FINN; Wiedinmyer et al., 2011), are used; FINN makes use of satellite and land coverage observations to estimate emissions from wildfires. Both EPA NEI-11 and FINN data are updated hourly to account for the diurnal cycle of their respective 125 emissions. Biogenic emissions of isoprene and monoterpenes are calculated online following Guenther et al. (1993). In-domain dust emissions are also calculated online following Zhao et al. (2010Zhao et al. ( , 2013b using the Goddard Global Ozone Chemistry Aerosol Radiation and Transport (GOCART) dust scheme from Ginoux et al. (2001). Short test simulations showed that surface dust concentrations were underpredicted. Therefore, the dust-tuning factor was increased from 1.0 to 1.2 in the six WRF-Chem experiments. 130

Non-dust aerosols
Most far-field cross boundary and initial condition chemistry for WRF-Chem is handled using the open-source software mozbc (NCAR, accessed 2018 at https://www.acom.ucar.edu/wrf-135 chem/download.shtml). This software uses simulated chemical output from the Goddard Earth Observing System version 5 (GEOS-5) model to generate the chemical initialization and lateral boundary condition files for WRF-Chem based on MOZART-4 (Model for OZone And Related Tracers version 4) simulations.
Chemical boundary tendencies are updated every 6 hours. Relevant to this study, mozbc is used to provide the lateral boundary and initial condition chemistry for many chemical constituents and aerosols, with the 140 exception of dust (see Section 2.3.2).
In this study, to avoid the issue of different dust bin cutoff sizes between MOZART-4 and MOSAIC, we ran a quasi-global (QG) WRF-Chem simulation to provide initial and cross-boundary dust in 145 our simulations. The QG WRF-Chem was ran at 1 o horizontal grid spacing with the identical options to the convective-permitting experiments except that convection was parameterized following Kain (2004). The

Experimental setup
Meteorological forcing at initialization and at the lateral boundaries of the convective-permitting experiments is provided from CFSR. Sea surface temperatures in the QG experiment and atmospheric 160 tendencies in the QG and convective-allowing experiments are updated every 6 hours. Convection is not parameterized in the 4-km experiments. A detailed listing of the physical parameterizations and physics packages used in the 4-km experiments are shown in Table 1.
The original control experiment was initialized on 26 September 2008 at 00:00 UTC and run through 1 August 2009. However, this experiment was found to underpredict snow water equivalent (SWE) 165 by several hundred millimeters due to underpredicted precipitation across the Rocky Mountains (not shown). A companion WRF simulation run without chemistry (NOCHEM) was found to significantly outperform WRF-Chem in simulating Rocky Mountain snowpack when compared to ground-based measurements; hence, a new set of WRF-Chem experiments was designed.
We restart our WRF-Chem simulations on 1 February 2009 00:00 UTC using surface energy and 170 hydrological fields from the NOCHEM restart file. Six "branch" WRF-Chem simulations are launched https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. from the new restart file to quantify the impacts of BCD across the Rocky Mountains (see Table 2). These experiments, to be run through 1 August 2009, consist of the following: 1.
CNT -the control experiment simulates both the SDEs and ARI associated with BCD. CNT 175 also includes the indirect effects associated with various aerosols, as the number concentration of activated aerosols is calculated based on the local aerosol characteristics in each grid cell.

2.
noSDE -this simulation is identical to CNT except that the deposition fluxes of BCD and snow BCD loading are set to zero in CLM/SNICAR. The deposition fluxes in the atmospheric component of WRF-Chem remain unchanged, effectively allowing BCD to vanish as they are 180 removed from the atmosphere. The ARI associated BCD remain in this perturbation experiment.

3.
noARI -this simulation is identical to CNT except that the mass contribution of BCD to the total aerosol optical properties is set to zero in the radiation code. Specifically, the contributions of BC, dust, and calcium to the atmospheric radiative effects are removed. 185

4.
noBCD -this simulation is identical to CNT except that the SDEs and ARI due to BCD are removed.

5.
noBCSDE -this simulation is identical to noSDE except that only the BC SDE is removed.
6. noBCARI -this simulation is identical to noARI except that only BC ARI are removed.
By examining the differences between the results of the six simulations, species-specific SDEs 190 and ARI associated with BCD can be quantified across 4 subregions of the Rocky Mountains shown in and (iv) the Southern Rockies; we consider elevations equal to or greater than 2,200, 2,400, 2,200, and 2,600 meters, respectively, across the 4 subregions. These regions were chosen because the water resources of these four areas depend heavily on the timing of local snowcover melt and orographic precipitation event 195 characteristics. These watershed-scale numerical estimates can serve to inform local policymakers, planners, and industries on the impacts of BCD effects. While BC SDEs and BC ARI are explicitly quantified by the difference between CNT and noBCSDE (noBCARI), it must be born in mind that dust SDEs (ARI) are taken to be the linear difference between noSDE and noBCSDE (noARI and noBCARI). https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License.  Daly et al., 1994).
PRISM fields are available on a ~4-km mesh. The spatial variability of monthly SWE is also evaluated in this study against data from the University of Arizona (UA). These data are generated with the methodology used to generate PRISM and are mapped to a ~4-km grid (Broxton et al., 2016). Model 210 simulated snow cover fraction (SCF) is evaluated against the High-resolution (0.05 o ) measurements from the MODerate resolution Imaging Spectroradiometer (MODIS) Aqua (Hall and Riggs, 2016).
Simulated BCD are compared to measurements from 23 Interagency Monitoring of PROtected Visual Environments (IMPROVE; Malm et al., 1994) network sites (see yellow triangles in Figure 1).
Simulated BC concentrations are approximated to be fine-mode elemental carbon (EC) from IMPROVE. 215 Although EC can be quite different from BC, this type of comparison has been made in other studies (e.g., Liu et al., 2012;. Simulated dust concentrations are approximated using the methods in Kevouras et al. (2007) by adding fine-mode soil to the difference between particulate matter (PM) having a dry-size diameter of less than 10 m (PM10) and PM having a dry-size diameter of less than 2.5 m (PM2.5). This approximation was found to be reliable at 9 inland rural IMPROVE sites in Malm et al. 220 (2007), at which the dust contribution to the coarse mode was 74-90%. The spatial distributions of February-through-July-averaged &( , #$% , and #&' for CNT and PRISM are shown in Figure 3, where large terrain-induced heterogeneity in temperature can be seen. For validation purposes, we exclude simulation results near CNT's lateral boundaries in an attempt to remove 240 gridpoints whose solutions are relaxed to coarse-resolution CFSR. All subsequent statistics make use of CNT data only within the black box shown in Figure  The spatial distributions of February-through-July-averaged precipitation rate, SWE, and SCF are shown in Figure 4. These hydrological fields do not correlate as well with observations as in the case of 2m temperature, with -values for CNT/PRISM precipitation rate, CNT/UA SWE, and CNT/MODIS SCF of +0.78, +0.90, and +0.82, respectively. Domain-averaged precipitation rate is overpredicted in CNT by 0.59 mm d -1 , with biases larger than 1 mm d -1 simulated locally at higher elevations. This wet bias is also 255 https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. evident in SWE and SCF comparisons, where CNT overpredicts SWE by more than 75 mm and SCF by 30% in several mountain chains. CNT simulates a mean wet bias that increases with increasing elevation; +, -values of +0.37, +0.20, and +0.38 are recorded for precipitation rate, SWE, and SCF, respectively.

BC and dust 260
The spatial distributions of in-domain February-through-July-averaged anthropogenic BC and dust surface emission rates are shown in Figures 5a and 5b, respectively. BC emissions are maximized within cities of the WUS such as Denver, Salt Lake City, and Las Vegas. Anthropogenic BC emissions also collocate with the federal interstate network. We note that, especially during summer months, a sizeable fraction of the BC emissions is associated with wildfires (not shown). Wildfires are scattered during our 265 simulation record (not shown), although they emit 100s to 1000s of milligrams per m 2 of BC into the atmosphere. Dust emissions on the other hand are maximized over the Mojave and Great Basin deserts as well as eastern Montana. As shown in Figure 5c and 5d, the largest burdens of BCD are simulated near or atop their respective emission sources. However, CNT emits over an order of magnitude more dust than BC, contributing to a large difference in February-through-July-averaged burdens (0.31 mg m -2 for BC, 9.9 270 mg m -2 for dust).

Radiative effects of BC and dust 300
Radiative effects of BCD are computed diagnostically following Ghan et al. (2012). For computation of the atmospheric radiative effects at the TOA and the surface, the radiation scheme is run interactively with all aerosol species factored into the calculation of beam transmittance (e.g. upwell shortwave flux at the TOA, downwelled longwave flux at the surface, etc.). The radiation scheme is then called iteratively in a methodology that sees, for example, BC excluded from the calculation of beam 305 transmittance variables. By subtracting the transmittance variables in the case where BC is excluded from the case including all aerosol species, the radiative effect of BC can be quantified. The same procedure is applied to quantify dust transmittance.

Atmospheric radiative effect 310
https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. May TOA RE of -1.05 W m -2 , but this region's TOA RE is comparable to other subregions during other months. The more negative TOA radiative effect due to dust across the Utah Mountains is due to the fact that this area sits immediately downstream of dust source regions. layer burdens of dust are ~2.5 orders of magnitude larger than those of BC (7.6 mg m -2 for dust compared to 0.0178 mg -2 for BC). ISREs are larger compared to CNT). Across the Utah Mountains and the Southern Rockies, the BCD ISRE peaks in April and May, respectively. The peak in ISRE during April across the Utah Mountains, occurring one month earlier than the other subregions, is due to maximized upstream dust emissions during April (not shown), which drives a fractionally larger dust-induced ISRE across this region. The reductions in SWE are driven by snow impurities that reduce snowpack albedo by as much as 0.02 across many elevated areas (Figure 8g). These albedo reductions are driven mainly by BC ( Figure 8h).

In-snow radiative effect
As the snowpack absorbs more heat, melting rates within the top snow are enhanced (not shown), leading to ice crystal growth of the underlying snow at the expense of liquid (Hadley and Kirchstetter, 2012). The geographic distribution of BCD SDE-induced anomalies can be tied to the mean normalized BCD burdens in the top snow layer (Figure 9). Burdens of BCD in the top snow represented in Figure 9 are divided by their respective means and are therefore unitless. The normalized in-snow burdens of BC are generally largest across western upslope regions of our subregions (Figure 9a). Different from in-snow BC 390 burdens, in-snow dust burdens generally decrease with northward extent across our domain (Figure 9b). This is due to the fact that southern mountain ranges sit directly downstream of dust emission sources, https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. while northern mountain chains do not. Additionally, the largest in-snow BCD mass mixing ratios are generally found across western upslope regimes where 1) aerosol deposition fluxes are maximized (not shown) and 2) snow burdens are relatively smaller to those further uphill (see Figure 4, middle row) 395 leading to maximized in-snow aerosol mass mixing ratios. The areas of maximum warming tend to be located within these moderately elevated upslope zones, while the largest SDE-induced reductions in SWE tend to occur at higher elevations downwind, where SWE is largest. The SWE results differ from those of , who found that the maximum reductions in SWE occurred across many basins of the WUS, especially across the northern subregions, although our result agrees with those of Qian et al. (2009). 400 It should be noted that WRF-Chem experiments undersimulate SWE and SCF across low-and moderatelyelevated regions such as the Snake River Basin (southern Idaho) and the Green River Basin (southwestern Wyoming). Simulated snow coverage may therefore have been too low for consequential SDE-induced perturbations to the meteorology across these areas, leading to the discrepancies between this study . 405

Timing of SDE
From Figure 10a it can be seen that BCD SDEs induce regionally averaged warming by no greater than 0.2 o C across our four subregions. The largest simulated warming occurs across the Utah Mountains in late April and Greater Idaho in mid-June (0.15 o C). BC almost universally heats the surface air (crosses, 410 Figure 10a), but we do note several instances where the dust SDE cools the surface air (hollow circles, Figure 10a). This is most probably the result of the assumption of linearity made in quantifying dust SDE as the difference between noBCSDE and noSDE (Section 2.3.2).
Peak SWE reductions are relatively well correlated with peak 2-m temperature increases and maximal ISRE values across the southern subregions. Peak SWE reductions of between 8 and 10 mm occur 415 in mid-April and mid-May across the Utah Mountains and Southern Rockies (between 4 and 5%), respectively. These reduction maxima, driven mainly by BC SDEs, occur concurrently with seasonally maximized ISRE values of +4 to +5 W m -2 ( Figure 11). It is noteworthy that SWE reductions across these southern subregions have a larger dust SDE-induced component than simulated across northern subregions.
As a percentage, negative SDE-induced SWE anomalies enlarge as the warm season progresses, with SWE 420 simulated reductions of ~65%, ~10%, ~7%, and ~50% occurring across Greater Idaho, the Northern Rockies, the Utah Mountains, and the Southern Rockies, respectively by mid-July. The progressive increase in SWE reduction percentages throughout the warm season is the result of decreased overall snowpack concurrent with increasing in-snow BCD mixing ratios as the snowpack ages, which increases the efficacy of BCD SDEs. 425 Northern subregions are not characterized by the same seasonal correlation between SWE, 2-m temperature, and ISRE. Maximized SWE reductions of 8 mm (4%) occur around 1 June across Greater Idaho, which occurs about a week after the ISRE maximum of (+2.9 W m -2 ) and ~3 weeks before the 2-m temperature maximum (+0.15 o C). This peculiarity is due to depressed snowmelt ( Figure 2b (Figure 10c), meaning that runoff anomalies are primarily driven by SDE-induced snowmelt. Here, runoff is defined to be the sum of surface and underground runoff; runoff from glaciers and lakes is neglected. Driven primarily by BC SDEs, runoff is mostly increased through late June across all four subregions. The largest increase in runoff occurs across the Northern Rockies (5.5 mm d -1 , a 90% change from CNT), which is characterized by the largest reductions in SWE. During July, negative 445 anomalies in runoff manifest, with the largest reductions simulated across the Northern Rockies (5.5 mm d -Southern Rockies in phase with precipitation increases across this subregion (Figure 10c). Although Qian et al. (2009) and emphasized results across basins, the dipole signature of runoff increases 450 followed by runoff decreases is consistent with our results, even though we primarily examine SDEs at higher elevations. SDE-induced perturbations to the weather and hydrological cycle generally stem from perturbations to the surface energy budget, affecting variables from the "bottom up." BCD ARI on the other hand seem to affect downwelled irradiance, atmospheric stability, and clouds in a way that affects the 455 underlying snow coverage at the surface, effectively altering variables from the "top down." Changes in clouds and low-level stability were minimal (less than 1% and 0.02 o C km -1 , respectively) and are not discussed.

Spatial patterns of ARI anomalies
From March through June, BCD ARI-induced 2-m temperature anomalies are minimal (less than 0.1 o C in magnitude) across the WUS (Figure 12a). However, BC and dust seem to impart anomalies of differing sign across the mountains, with cooling (warming) simulated across the higher elevations due to BC (dust). This general spatial pattern is correlated with SWE increases due to BC (Figure 12e) and SWE 465 decreases due to dust (Figure 12f). The cooling and warming patterns associated with BC and dust, respectively, can be tied to deficits and surpluses in the surface energy budget imparted by these aerosols' atmospheric radiative effects (to be discussed in section 5.3). Figure 12d shows that BC ARI-induced SWE increases tend to exceed reductions in SWE associated with dust ARI. These overall ARI-induced increases in SWE compete against, but do not 470 exceed, BCD SDE-induced reductions in SWE (Figure 8d). Total simulated precipitation (rain+snow) anomalies due to BCD ARI are generally less than 0.2 mm d -1 , while those due solely to snow are less than 0.1 mm d -1 (not shown). Additionally, no precipitation variable shows any discernable weekly, monthly, or seasonal trend. This finding lends further credence to the idea that simulated ARI-induced changes to SWE manifest from changes to the surface energy budget through BCD ARI. 475 https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. Figure 12g shows that BCD ARI imparts a negative surface RE. Curiously, BC and dust have surface REs that are opposite in sign (Figures 12h,i). The opposing sign is attributable to the differing microphysical properties of BCD aerosls (discussed in Section 5.3). Figure 13a shows the seasonality of competing 2-m temperature anomalies with BC (dust) ARI inducing surface cooling (warming). Maximum BCD ARI-induced temperature perturbations tend to occur from late May into early June. Interestingly, the strongest dust ARI-induced warming (+0.1 o C in late May) tends to occur across the Northern Rockies, even though this region lies further away from dust emission sources. As will be discussed, this region tends to have relatively lower burdens of super-micron dust 485 particles compared to the other subregions. As larger (smaller) dust particles tend to absorb (scatter) incoming and upwelling sunlight, this region tends to see more downward-scattered insolation due to the larger fraction of sub-micron dust particles, which results in less surface dimming. BC on the other hand absorbs in the atmosphere, leading to surface dimming and weak surface cooling. Otherwise, the strongest combined (BC+dust) surface cooling is generally simulated to occur across the southern subregions prior to 490 mid-May, but all subregions tend to see BCD ARI-induced cooling throughout the simulation period.

Timing of ARI 480
Overall cooling by BCD ARI is accompanied by increased SWE amounts across the four subregions ( Figure 13b). Driven by BC ARI, SWE increases of 5 mm (2%) are simulated across the Northern Rockies during mid-May. By percentage however, the largest SWE increases of ~3 mm (~4%) are simulated during the first week of May across the Utah Mountains. Across the southern subregions, the 495 ARI-induced increases in SWE occur earlier in the spring relative to the northern subregions. Meanwhile, the percentage increases in SWE across all four subregions generally increase with time well into the summer, with increases in excess of 10% across the Utah Mountains, the Southern Rockies, and the Northern Rockies by late July (not shown).
Overall, the increase in SWE through ARI is predominately due to BC; however dust ARI does 500 contribute to SWE increases on the order of 1 mm during mid-May and early June across the Utah Mountains and the Southern Rockies, respectively. Interestingly, the Northern Rockies see competition between BC and dust ARI in modulating SWE anomalies. Specifically, BC (dust) increases (decreases) https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. SWE by 8 mm (4 mm). Differences in these aerosols' optical properties lead to differences in how these aerosols perturb the surface energy budget while they reside in the atmosphere. These differences are 505 discussed in Section 5.3; however, it is noteworthy that the changes to snowpack are not due to simulated ARI-induced changes in precipitation (not shown). Specifically, ARI-induced precipitation anomalies are similar to those induced by SDE in that they show no discernable trend and generally do not exceed 0.3 mm d -1 in magnitude. Because precipitation anomalies are on the order of tenths of millimeters per day, and SWE anomalies may be more than 10 mm locally, this indicates that ARI-induced changes in SWE are due 510 mostly to changes in how the snowpack melts with time rather than changes to precipitation.
With SWE increased from April onward due to BCD ARI, runoff tends to decrease across all four subregions prior to mid-May (Figure 13c). Following 1 June, runoff anomalies become less negative and even positive across the four subregions, a pattern opposite to that of BCD SDE (Figure 10d). BC ARI tend to drive a majority of the runoff decreases prior to 1 June and promote increased runoff deeper into the 515 summer. Dust ARI on the other hand has the opposite effect on runoff to that of BC ARI, increasing runoff through mid-May and decreasing runoff after 1 June.

Main differences between SDE and ARI in affecting snowmelt
Overall, BC SDE incites anomalies opposite to those of BC ARI. The former heats the snow from 520 within, introduces a positive radiative perturbation at the surface, and increases melting through surface albedo reductions. The latter effect preserves the snowpack through aerosol dimming at the surface.
Specifically, BC absorbs direct incoming sunlight, leading to warming of the TOA (Table 3). The absorption of light by BC within the atmosphere reduces the amount of insolation at the earth's surface, resulting in a surface RE associated with BC of between -0.5 W m -2 and -1 W m -2 (see Figure 12h and 525 Figure 13d; cross symbols). Furthermore, this dimming appears to be maximized across the southern subregions during mid-May.
Changes to the snowpack brought about by dust ARI and dust SDE are generally smaller than those brought about by BC SDE and BC ARI. However, a notable instance in which dust ARI brings forth SWE reductions of 4 mm around 1 June across the Northern Rockies (Figure 13b) can be explained by 530 considering the microphysical characteristics of dust particles across the subregion. Regional averages of https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License.
dust-(and BC-) induced anomalies were computed over areas characterized by high albedos. Smaller dust particles are generally scattering in nature compared to larger dust particles in shortwave bands (Tegen and Lacis, 2012). In a single-layer radiative transfer model of the atmosphere, the downward shortwave solar flux ↓ at the Earth's surface can be expressed as: 535 where 0 is the solar constant at a particular time of day at a latitude characterizing the Northern Rockies and subject to molecular scattering, 5 is the fraction of 0 scattered by aerosols in the atmosphere, and & is the fraction of 0 absorbed by aerosols in the atmosphere. The upwelled solar flux immediately above the earth's surface ↑ can thus be found by multiplying (1) by the surface albedo : This upwelled radiation will be scattered back and forth between the scattering aerosol layer and the 540 surface in an iterative process proportional to 5 and : (3) Here, n is the number of times the incident beam is scattered between the surface and the aerosol scattering layer. <=<&> ↓ − 0 is negative. However, by (3), it is clear that in a case where & is negligible, <=<&> ↓ − 0 becomes less negative as approaches unity. For atmospheric dust particles residing over the high-albedo surface of the Northern Rockies, this means that there will be a higher chance of shortwave absorption at 545 the surface. The reduced negativity of the dust surface shortwave RE, together with the dust longwave warming (Figure 12i and Figure 13d), explains why dust aerosols contribute to snowpack reductions across the Northern Rockies. The physical process described here is similar to that noted in Stone et al. (2008) who examined the atmospheric REs of wildfire smoke across northern Alaska's high albedo surface. 550

Consequences of increased dust
The mean [dust] was undersimulated in all WRF-Chem experiments by 63% compared to IMPROVE (Section 3.2). Even in the DTF=2 experiment, [dust] was still underpredicted by 43%. Figure   14 shows the consequences of increasing dust emissions by 60% from March through June. ISRE is increased by 0.2 W m -2 in many areas, especially across the southern subregions. Additionally, the northern 555 https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. subregions see ISRE increases of 0.05 to 0.1 W m -2 (Figure 14a). Across Greater Idaho and the Northern Rockies, this represents a doubling of the ISRE in the DTF=2 experiment compared to CNT.
As discussed in the previous section, in-atmosphere dust can warm the surface through shortwave and longwave ARIs over high-albedo surfaces. Figure 14b shows the difference in the net surface RE between the DTF=2 experiment and CNT. A positive differnence of between 0.01 and 0.05 W m -2 is 560 prevalent across the four subregions, especially across the Northern Rockies and Southern Rockies. Combined with BC ARI dimming, the dust ARI offsets snow darkening reductions in SWE, resulting in relatively smaller net perturbations to the surface water budget.

Conclusions 595
Using seven "branch" WRF-Chem experiments with a horizontal resolution of 4 km, the SDEs of BCD were quantified across four subregions of the WUS for water year 2009. These aerosols' ARI were also examined and quantified.
It was found that BC surface concentrations in the WRF-Chem control experiment (CNT) were well simulated compared to observations from IMPROVE, while dust surface concentrations [dust] were 600 undersimulated by 63%. An additional simulation was run with domain-resolved dust emissions increased by 60% (DTF=2), but this simulation still underpredicted [dust] by a factor of 43% compared to IMPROVE. It was found that CNT generally overpredicted precipitation compared to PRISM, snow cover compared to MODIS, and SWE compared to UA at all elevations; the wet bias increased with increasing elevation. Meanwhile, CNT simulated a warm (cold) bias at lower (higher) elevations. The biases/elevation 605 Pearson values exceeded +0.3 for hydrologic reference variables and fell below -0.65 for the temperature variables.
It was found that BCD SDEs generally reduce springtime and summertime SWE while ARIs generally increase SWE during this same time period. However, SWE reductions due to the SDE, typically on the order of 2% to 5%, far exceeded the increases induced by ARIs. SWE changes by these BCD effects 610 led to changes in runoff. Specifically, SDEs bring forth runoff increases of 1-5 mm d -1 before July and runoff reductions of a similar magnitude during July. Runoff increases and decreases were largest across the Northern Rockies, with runoff increases of 95% during springtime and early summer preceding gradual decreases of 5% through the first half of July. BC ARI on the other hand dims incoming solar radiation by 0.5-1.0 W m -2 , depressing snowmelt in the late spring. The resulting ARI-induced perturbation to runoff is 615 generally opposite to that associated with SDE. To emphasize, the changes to surface hydrology across the WUS are driven primarily by SDE, not ARI, and simulated SWE anomalies brought forth by BCD were not due to precipitation anomalies.
It was found that BC ARI dims the surface, resulting in a negative surface RE. Depending on the subregion however, dust ARI could actually incite a positive surface RE. We believe this to be the result of 620 differing optical properties associated with dust aerosols of differing size. Across subregions in relatively close proximity to dust emissions sources, larger, more absorbing dust particles dim the surface, preserving snowpack through a negative surface RE. Across regions further away from dust emission zones, namely the Northern Rockies, dust particles tend to be smaller and more scattering. Across this high-albedo surface, downwelled solar radiation is reflected upward and subsequently backscattered to the surface by 625 the scattering dust layer in an iterative process. Suppressed dimming of incoming solar energy, in tandem with a positive longwave RE due to dust ARI, warmed the surface and reduce SWE.
Indeed, BC-induced effects on WUS meteorology and hydrology for water year 2009 were generally simulated to be larger than those induced by dust. We conclude however by noting that there were non-negligible changes in SWE in the instance that dust emissions were increased. Specifically, SWE 630 reduction enhancements of several millimeters were simulated. While our results were consistent with Wu et al. (2018), we note that observed [dust] was estimated following Kevouras et al., (2007) and this estimation is subject to uncertainties therein. There remains a possibility that dust effects may be quite comparable to those of BC, although further simulations are required to answer this question. More generally, BCD SDEs and ARI can impart significant perturbations of WUS weather and hydrology and 635 corroborate the results of coarser resolution GCMs. Future studies should focus on increasing the domain size in order to quantify larger synoptic-scale circulation changes associated with BCD effects; this coupling has been noted in previous GCM experiments (e.g. Rahimi et al. 2019), but the domain size used https://doi.org/10.5194/acp-2019-998 Preprint. Discussion started: 9 January 2020 c Author(s) 2020. CC BY 4.0 License. in this study was too small to resolve these potentially significant larger-scale responses. Additionally, because convection was not parameterized, an opportunity exists to quantify how orographically forced 640 precipitation events are impacted by BCD SDE. Specifically, this output data can be used to quantify how updraft speeds, convective updraft area, and storm energetics change as a function of BCD effects, which will be topics of future study.  , 9, 1959-1976, doi: 10.5194/gmd-9-1959-2016, 2016. 860         c. d.

ARI-induced perturbations by BCD
a. b. c.