Global Radiative Impacts of Mineral Dust Perturbations Through Stratiform Clouds

Airborne mineral dust influences cloud occurrence and optical properties, which may provide a pathway for recent and future changes in dust concentration to alter the temperature at Earth's surface. However, despite prior suggestions that dust‐cloud interactions are an important control on the Earth's radiation balance, we find global mean cloud radiative effects to be insensitive to widespread dust changes. Here we simulate uniformly applied shifts in dust amount in a present‐day atmosphere using a version of the CAM5 atmosphere model (within CESM v1.2.2) modified to incorporate laboratory‐based ice nucleation parameterizations in stratiform clouds. Increasing and decreasing dustiness from current levels to paleoclimate extremes caused effective radiative forcings through clouds of +0.02 ± 0.01 and −0.05 ± 0.02 W/m2, respectively, with ranges of −0.26 to +0.13 W/m2 and −0.21 to +0.39 W/m2 from sensitivity tests. Our simulations suggest that these forcings are limited by several factors. Longwave and shortwave impacts largely cancel, particularly in mixed‐phase clouds, while in warm and cirrus clouds opposite responses between regions further reduce each global forcing. Additionally, changes in dustiness cause opposite forcings through aerosol indirect effects in mixed‐phase clouds as in cirrus, while in warm clouds indirect effects are weak at nearly all locations. Nevertheless, regional forcings and global impacts on longwave and shortwave radiation were found to be nonnegligible, suggesting that cloud‐mediated dust effects have significance in simulations of present and future climate.


Introduction
Airborne mineral dust influences Earth's radiative balance through a variety of mechanisms involving clouds, complementing the direct scattering and absorption of radiation by dust particles (dust direct effects). Clouds both scatter incoming solar radiation and block outgoing longwave radiation, with mineral dust impacts on cloud occurrence and optical thickness influencing the balance between these effects. Mineral dust affects cloud presence and properties through thermodynamic responses to dust layer heating (dust semidirect effects) and dust's roles in cloud droplet and ice crystal formation (dust indirect effects). Due to the physical and radiative properties of mineral dust and its abundance as one of the atmosphere's most globally prevalent aerosol in terms of mass (Textor et al., 2006), it is responsible for a significant proportion of aerosol scattering and absorption of solar radiation (Heald et al., 2014;Samset et al., 2018) and is likely the aerosol species most commonly active as ice nucleating particles (INPs) for ice crystal formation in both cirrus and mixed-phase clouds (defined below) (Cziczo et al., 2013;Murray et al., 2012).
Aerosol indirect effects present some of the largest uncertainties in simulations of present and future climate (Boucher et al., 2013) due to their complex nature, which includes different pathways in warm clouds containing liquid water droplets and cold clouds consisting either purely of ice crystals (cirrus) or both ice crystals and liquid droplets (mixed-phase clouds). Through its role as INPs, mineral dust is directly involved in heterogeneous ice nucleation, the process by which ice crystals form on the surface of particles, as compared to homogeneous nucleation, which involves no such surface (Vali et al., 2015). Heterogeneous nucleation is traditionally divided between nucleation on INPs immersed in a droplet (immersion freezing), on INPs where freezing occurs as vapor condenses on the surface (condensation freezing), at interfaces where INP contacts droplets (contact freezing), or by deposition of vapor on the INP (deposition nucleation) (Vali et al., 2015). Immersion and condensation freezing events are difficult to distinguish in observations, so models tend to combine these mechanisms. In mixed-phase clouds, immersion freezing is thought to be the dominant nucleation mechanism (de Boer et al., 2011;Hoose et al., 2010), while homogeneous nucleation does not occur because temperatures are too high. In cirrus clouds, heterogeneous nucleation was traditionally attributed to deposition nucleation, but this is now thought to be initiated by the freezing of liquid water trapped within pores in aerosol surfaces (the pore condensation and freezing [PCF] mechanism) (David et al., 2019Marcolli, 2014). Analyses of satellite retrievals suggest that dust acting as INPs influences ice crystal number in mixed-phase (Zhang et al., 2018) and cirrus clouds  and enables glaciation of mixed-phase clouds (Tan et al., 2014). Heterogeneous ice nucleation enabled by the presence of dust in cirrus clouds can create a sink of water vapor that prevents homogeneous nucleation from occurring to such an extent that dust may be either positively or negatively correlated with ice crystal number density and size depending on conditions (Kuebbeler et al., 2014;Zhao et al., 2018;Zhou et al., 2016), with regional studies of the influence of dust being potentially unrepresentative of the global effect.
Global mineral dust concentrations have varied immensely throughout time. Ice cores have evidenced that during the Last Glacial Maximum (LGM, 21,000 years ago) global mineral dust concentrations were two to five times higher than during present-day (Kohfeld & Harrison, 2001), while global dust concentrations during the mid-Pliocene warm period concentrations are unknown but thought to be far lower than current levels (Sagoo & Storelvmo, 2017). Estimates of more recent dustiness suggest that global concentrations have doubled during the Anthropocene (Hooper & Marx, 2018) and possibly over the 20 th century alone (Mahowald et al., 2010). Projections of future dustiness are difficult to make due to dependence on changes in surface winds (Evan et al., 2016;Kok et al., 2018;Pu & Ginoux, 2018), precipitation (Kok et al., 2018;Pu & Ginoux, 2018), surface bareness (Pu & Ginoux, 2018), and land use in arid regions (Ward et al., 2014). CMIP5 projections of dust burden changes over the 21st century have ranged from −17% to +23%, with a +5.2% model mean, under a Representative Concentration Pathway (RCP) 8.5 scenario (Allen et al., 2016). An understanding of the global impacts of past or hypothetical dust perturbations on clouds and radiation may guide research on the extent to which such changes can potentially intensify or offset anthropogenic climate change.
In this study we simulate large, uniformly applied mineral dust perturbations with the Community Earth System Model (CESM) version 1.2.2 and assess the radiative impacts involving dust-cloud interactions. As recent advances have occurred in the understanding of indirect effects in cirrus and mixed-phase clouds, the default CESM does not represent the best available knowledge of microphysical responses to dust changes. To address this limitation, we have modified the model's atmospheric component to include contemporary laboratory-based parameterizations of dust ice nucleating potential from Ullrich et al. (2017) (hereafter referred to as Ull17). We assess the influence of dust changes on cloud properties and the ensuing radiative forcings, which we separate into four mechanisms using an original decomposition method: cirrus indirect, mixed-phase cloud indirect, warm cloud indirect, and semidirect. Sensitivity tests were conducted to test a range of plausible scenarios that could affect the strength of these forcings. The result is a comprehensive bound on the radiative forcings that would be induced by the simulated dust perturbations and incorporates the most current knowledge of aerosol-cloud interactions.

Core CESM Simulations
We ran 5-year climate model simulations with varied dust emissions in sets having dust-cloud and dust-climate linkages variously enabled and disabled, and compared simulations to quantify the radiative impacts of dust emission differences attributable to each discussed mechanism. Simulations are listed in Table 1, with analysis methods described in section 2.4. All simulations used CESM v1.2.2 with the Community Atmosphere Model v5 (CAM5) (Hurrell et al., 2013) at 1.9°× 2.5°resolution, with the threemode Modal Aerosol Module (MAM3)  and Morrison and Gettelman (2008) stratiform microphysics (hereafter MG08). We modified MG08 to make ice nucleation rates in cirrus and mixedphase clouds respond to mineral dust concentrations in tune with laboratory observations. Simulated dust perturbations were extreme and applied uniformly (as set fractional changes at all months and locations globally), sea surface temperatures were fixed, and winds were nudged toward the NASA Modern-Era Retrospective Analysis for Research and Applications (MERRA) reanalysis (Reinecker et al., 2011) over the period 2004-2008 using a 6-h relaxation time scale. This setup enabled clean signals of dust-cloud forcings with moderate computational cost but rendered us unable to account for climate feedbacks involving atmospheric and oceanic circulations, including the potential for impacts on winds at surface and aloft to feed back on dust emissions and transport. While temperature and other meteorological variables were available in the reanalysis, due to potential biases (Zhang et al., 2014) we chose only to nudge winds. All simulations had a 1-month spin up (December 2003), not used in calculations, in order to initiate the dust cycle, and used corrected dust optical and hygroscopic properties (Albani et al., 2014).
Mineral dust emissions were prescribed as an annual cycle from emission files using the methodology of Lou et al. (2017). Dust emission files for present day were created from output of an initial 13-month CESM simulation (with the first month discarded), using the default dust model (Zender et al., 2003) with modified soil erodability and dust size distribution (Albani et al., 2014). For this simulation, winds were nudged toward MERRA reanalysis for the year 2004, and dust emission strength was tuned such that for a full nudged run from 2004-2008 the global average dust aerosol optical depth (AOD) was within 0.030 ± 0.005, as has been estimated from satellite observations of that period (Ridley et al., 2016). In our present-day run dust AOD averaged to 0.033 globally over the full 5-year period. For simulations for low dust and high dust cases, the same emissions were fed in with all values multiplied by 1/10th and 10/3rds, respectively. The dust levels used here are similar to those assumed for the mid-Pliocene warm period and LGM in Sagoo and Storelvmo (2017), respectively. While these represent far more extreme cases than those projected for near-future climate, their use reflects the atmosphere's sensitivity while sufficiently perturbing the model for a clear signal to develop. We depend on the model to appropriately simulate the vertical distribution of dust, as we have only constrained the columnar dust AOD rather than proportions at individual levels.
A number of changes were made to the model's stratiform microphysics routines to better estimate cirrus and mixed-phase indirect effects. These are described in section 2.2. The simulated ability of mineral dust Note. For each set three simulations were run: present-day dust, low dust, and high dust. "Def" implies the default or most standard setting.

10.1029/2019JD031807
Journal of Geophysical Research: Atmospheres to act as cloud condensation nuclei (CCN) in warm clouds is kept as in the default model, using a scheme (Abdul-Razzak & Ghan, 2000) based on Köhler theory. In this scheme each aerosol mode is considered internally mixed. Aerosol activation rates are calculated based on mass-weighted mean hygroscopicity and are influenced by the number and size of aerosol in each mode as well as by vertical velocity. Dust radiative impacts through direct scattering and absorption are represented in our simulations, which include updated dust optical properties (Albani et al., 2014), but are not the focus of this study and are covered elsewhere (e.g., Kok et al., 2017). Other simulated dust radiative effects involve dust impacts on water vapor concentrations and land surface temperature as well as surface albedo effects of dust deposition on snow. We hereafter refer to these and any other dust radiative effects distinct from dust-cloud and direct dust effects as remaining dust effects.

Modifications of Ice Nucleation in CESM
We replaced all ice nucleation schemes affecting stratiform cold clouds with the newer empirical parameterizations of Ull17. The Ull17 parameterizations are based on 11 years of measurements from the AIDA cloud chamber and estimate ice crystal number for both deposition nucleation on dust in cirrus clouds and immersion freezing on dust in mixed-phase clouds. The parameterizations calculate activated ice crystal number given meteorological conditions and dust surface area and abundance, building on the concept of ice nucleation active surface (INAS) density (Connolly et al., 2009;Hoose & Möhler, 2012). Separate parameterizations in the Ull17 study pertaining to ice nucleation on soot-of highly uncertain importance (Kanji et al., 2017) -were not used. We treated dust INPs as globally uniform in composition (as in the parameterization) and, in most simulations, as freshly emitted rather than partially coated by other aerosol species. This second assumption could overestimate the ice nucleating ability of aged dust in the cirrus temperature range . Dust's impacts on ice nucleation in convective clouds are less well studied than in stratiform clouds and are not directly represented in our simulations.
For ice nucleation in stratiform mixed-phase clouds, we replaced all previously used schemes (Bigg, 1953;Meyers et al., 1992;Young, 1974) with the Ull17 dust immersion scheme. The Ull17 parameterization is solely a function of temperature and mineral dust, as water vapor in mixed-phase clouds is assumed to always be at saturation with respect to liquid water. We incorporated Ull17 into MG08 as a source of newly activated ice crystals enabled by the adiabatic cooling of dust-laden air parcels within sub-grid scale updrafts. Our implementation first calculates the adiabatic cooling of sub-grid scale updrafts during each time step given an updraft speed empirically derived from the model's estimate of turbulent kinetic energy (TKE). Conditions before and after this cooling are then fed into the parameterization with the difference in calculated density of heterogeneously frozen crystals deciding the ice crystal source. By comparison, the Meyers scheme, which is the most dominant local ice source in default MG08, is aerosol independent and is implemented so that ice crystal number is pushed toward a number calculated for the given temperature if the current value is below the calculated one, a setup that (a) prevents nucleation from occurring in situations when a large amount of crystals from other sources (crystals sedimenting from cirrus clouds or convective detrainment) make their way into stratiform mixed-phase clouds and (b) enables fast nucleation in conditions of strong ice crystal sinks. Whereas Meyers accounts for deposition nucleation and assumes a 10-micron diameter of new ice crystals, Ull17 assumes immersion freezing to be the dominant nucleation process in mixed-phase clouds; hence, our implementation of Ull17 takes water involved in stratiform mixed-phase cloud ice nucleation from the liquid phase and sets the size of new crystals to match the mass of co-located liquid droplets. Previous implementations of aerosol-sensitive immersion freezing parameterizations in CESM (Sagoo & Storelvmo, 2017;Tan & Storelvmo, 2016) had been similar to the Meyers scheme in that they also pushed ice number toward a calculated density and took mass for new ice crystals from the vapor phase without limitations from droplet availability. Note that we assume all local nucleation in mixed-phase clouds to come from Ull17 and do not include parameterizations for deposition nucleation or contact freezing in mixed-phase clouds when using this parameterization. To match our changes to heterogeneous ice nucleation in mixed-phase clouds, we also strengthened homogeneous freezing nucleation to freeze all cloud droplets that adiabatically cool to −40°C within a time step rather than only those within grid cells already at −40°C as in default MG08. Our simulations did not take into account influences of coatings on dust when calculating mixed-phase cloud ice nucleation rates. Such coatings have not been found to cause great influence in this cloud regime (Kanji et al., 2019).
Due to the overabundance of cloud ice relative to cloud liquid at mixed-phase temperatures (roughly −38°C to −5°C), the process of ice crystals growing at the expense of liquid droplets (the Wegener-Bergeron-Findeisen, or WBF, process) was weakened to 1% of its initial value as in Sagoo and Storelvmo (2017). We represented this in the model as cloud liquid and ice being truly in contact only in 1% of clouds that contain both liquid and ice, with the rest being pockets of ice only and liquid only. In the remaining cloud ice, existing ice crystals were enabled to grow at the expense of water vapor unlike in default MG08, to avoid sudden jumps in ice depositional growth rate from zero when the last droplet is glaciated. To further account for an overabundance of cloud ice phase compared to liquid, we modified the thermodynamic phase of convective detrainment. Ice nucleation in CAM5 deep and shallow convection is highly simplified, with the phase ratio of new condensate (cloud liquid to cloud ice) being linearly scaled to temperature and independent of aerosol concentration. Here we limited the detrainment of ice to temperatures colder than −35°C, as described in Tan and Storelvmo (2016). Figure 1a shows that cloud phase in our simulation of present-day climate roughly matched that of satellite observations along relevant isotherms, while Figure 1b shows a small to moderate underestimate of ice crystal number in this simulation.
For ice nucleation in cirrus clouds, the parameterizations of the default model scheme (Liu & Penner, 2005) were replaced by Barahona and Nenes (2009) (hereafter BN09), and the heterogeneous nucleation parameterization in BN09 was replaced with the deposition parameterization of Ull17. For a smaller group of runs the PCF-based parameterization of David et al. (2019) (hereafter Dav19) was instead used as the heterogeneous parameterization. The Ull17 parameterization estimates the amount of ice crystals as a function of temperature, humidity, and mineral dust surface area and number density. The Dav19 parameterization is not directly temperature dependent, only being a function of humidity and dust number, with separate functions for fine and coarse mode dust. Both parameterizations were embedded within the BN09 ice nucleation scheme, which simulates the competition between heterogeneous nucleation and the homogeneous nucleation of sulfate droplets. BN09 calculates an expected number of ice crystals given temperature, updraft speed, saturation level, and INP and sulfate amounts and is integrated with MG08 to add ice crystals where less than calculated are present. In most simulations the mass transfer coefficient of water vapor onto ice (deposition coefficient) was kept at the relatively high value of 1.0 (Zhou et al., 2016) to approximately match ice crystal number with observations ( Figure 1b). The sub-grid scale vertical velocities fed into the cirrus nucleation scheme were a function of TKE, as in several studies (e.g., Gettelman et al., 2010;Lohmann et al., 1999). We treated vertical velocity as a normal distribution around the resolved vertical wind with the standard deviation being a constant multiplied by the square root of TKE. Our choice of 0.6 for this constant was found to roughly match global estimates of vertical velocities based on models and observations The green solid line is from a simulation with present-day dust levels while blue and red represent low and high dust scenarios, respectively. The black lines represent satellite retrievals, with SLF data from CALIOP as in Figure S1 of Tan and Storelvmo (2016) in Panel (a) and ice crystal number data from Sourdeval et al. (2018) in Panel (b). Solid lines represent global averages while dashed lines in Panel (b) denote spatial medians with shaded areas showing the 25th and 75th percentiles of the present-day simulation and observations. All values in Panel (b) are in-cloud measures that incorporate only fully glaciated clouds, as retrieval of ice crystal number is difficult in clouds containing liquid. Model output shown is from simulation set C1 as listed in Table 1.

10.1029/2019JD031807
Journal of Geophysical Research: Atmospheres (Barahona et al., 2017), albeit with higher vertical velocity variability in the tropics and lower elsewhere. We then divided the distribution into six parts and set the BN09 scheme to estimate the integration over the distribution via Gauss-Legendre quadrature. In mixed-phase clouds, we based the updrafts on the same vertical velocity distribution used in cirrus clouds but fed in a single average updraft in each time step rather than the distribution itself. This average updraft velocity was set to the standard deviation multiplied by √(2/π) (Fountoukis et al., 2007), and we assumed it to be operating over half of the cloud area with the remainder undergoing downdrafts and no new ice nucleation. Simulations neglected the potential weakening of cirrus nucleation on dust due to coatings  except in specified sensitivity test runs. All simulations neglected the potential for cloud processing of INPs to strengthen cirrus ice nucleation Umo et al., 2019;Wagner et al., 2016), as our model did not keep track of which aerosol had already been used for ice nucleation and the effect is not included in most current parameterizations.

Sensitivity Test Simulations
Several uncertainties exist in simulations of dust-cloud-climate interactions. Two major factors are that (a) parameterizations of heterogeneous ice nucleation activity on dust do not fully agree with one another and (b) climate model constants and input variables related to dust radiative properties and nucleating ability are often highly uncertain and poorly represented. To account for these uncertainties we performed a number of sensitivity simulations and recalculated the associated dust radiative effects.
Parameterization disagreements may represent differences in aerosol samples, as little is known about how representative any laboratory aerosol is to those in the atmosphere (which may vary by source and be aged by other species) or differences in equipment and experimental design. We reran tests of dust indirect effects replacing the Ull17 immersion nucleation parameterization with that of DeMott et al. (2015) (DM15 hereafter) and, separately, replacing the Ull17 deposition parameterization with that of David et al. (2019) (Dav19), both with the same implementations as described. DM15 calculates activated ice crystal number in mixed-phase clouds based on field observations using a continuous flow diffusion chamber (CFDC) (Rogers et al., 2001), with the number of dust particles having >0.5-μm diameter as input. Dav19 estimates heterogeneous ice nucleation in cirrus clouds as a function of total dust number concentration based on measurements from the Zurich Ice Nucleation Chamber (ZINC) CFDC instrument (Stetzer et al., 2008) and theory under the basis that heterogeneous ice nucleation at cirrus temperatures acts primarily by the PCF mechanism. We additionally reran simulations using the Ull17 deposition scheme modified  to have dust act as weaker INPs to emulate the potential impacts of dust being atmospherically aged uniformly.
To address input uncertainties, we performed simulations with the default Ull17 schemes while halving and doubling the vertical velocities fed into each microphysical scheme (droplet activation and ice nucleation in both mixed-phase and cirrus clouds), and an additional set of simulations with the deposition coefficient set to 0.1 (at the lower end of the observational range, Zhou et al., 2016) and cirrus updrafts weakened to offset the impact on ice crystal number. Semidirect effects are susceptible to uncertainty in dust's ability to absorb solar radiation, so we test the range of imaginary refractive index values from Kok et al. (2017).

Partitioning of Cloud Forcings
In all simulations, the shortwave and longwave radiative balance in each model column was calculated with and without clouds, with the differences being the cloud forcing (Ramanathan et al., 1989). To distinguish the indirect and semidirect impacts of mineral dust, we ran simulations both with and without mineral dust scattering and absorption effects on. We also added clean-air diagnostics (calculated as if aerosols were not present), as in past studies of anthropogenic aerosol effects (Ghan, 2013;Ghan et al., 2012) though with diagnostics added to account for both shortwave and longwave effects rather than shortwave only. This allowed us to separate aerosol direct effect forcings from forcings due to other factors independent of clouds (primarily surface albedo effects via deposition onto snow and ice, water vapor impacts, and differences due to altered land surface temperatures).
To partition the forcings due to varied dust levels into that attributable to each studied mechanisms, we ran sets of simulations having dust direct effects and dust-cloud linkages variously enabled and disabled. Shown in Table 1 are the core simulation sets ("C") necessary for our results and additional simulations for sensitivity tests ("T"). In Equation 1 we present the subtractions between simulations used to calculate each associated effective radiative forcing (ERF) for our core estimate calculations, noting simulation sets in subscripts (e.g., "C1"). These calculations were done for both shortwave and longwave top-of-atmosphere (TOA) forcings, with net forcings calculated as their sum. Note that cloud forcings are calculated by subtracting radiative diagnostics neglecting clouds (marked clear) from the actual radiative forcings, while dust ERFs are calculated as differences between simulations having different dust levels ("Δ", representing for instance "high dust"-"present day dust"). To calculate the semidirect effect we used diagnostics neglecting all aerosol (marked clean), in order to remove an artifact of the direct effect forcing, and hence, the sums of these cloud forcings do not equal the total cloud forcings from the model. We further removed a second, smaller forcing from the semidirect effect, representing dust impacts on the cloud forcing difference with and without aerosol effects in simulations without dust direct effects, as this also did not appear to be directly relevant. A brief description of these "remaining" effects is given in supporting information Text S2 and shown in Figures S1 and S2. For simulations with dust direct effects off (C2-4 in the core simulations), we disabled all scattering and absorption by dust along with the radiative effects of redistributed aerosol water onto dust. For those with heterogeneous dust nucleation in cirrus off (e.g., C4), all cirrus nucleation was homogeneous nucleation of sulfate. For simulations where dust nucleation was off in mixed-phase clouds (e.g., C3 and C4), we instead turned on the Meyers deposition nucleation scheme, which is aerosol independent, to prevent ice crystal numbers from reaching unrealistically low values. Each sensitivity test was conducted by replacing one core simulation set in Equation 1 with a related sensitivity set: T1 and T2 replace C4 (warm cloud indirect effect tests), T8-10 replace C2 (mixed-phase cloud indirect tests), T3-7 replace C3 (cirrus indirect effect tests), and T11 and T12 replace C1 (semidirect effect tests).
The utilized method of isolating dust-cloud effects by removing dust effects one after another is not perfect, as the removal of each effect slightly alters the base state of clouds involved in the calculation of unassociated forcing mechanisms. Hence, the removal of dust direct effects for the calculation of indirect effects, indirect effects in mixed-phase clouds for their calculation in cirrus and warm clouds, and indirect effects in cirrus for their calculation in warm clouds may induce biases that cause total dust-cloud effects to be partly attributed to the wrong mechanisms. In Figures S6-S12 we show that as dust-cloud effects were removed the base states of evaluated clouds remained similar, letting us conclude that calculations were not markedly compromised by this issue. Figures S6-S11 show cloud microphysical variables in simulations with present-day dust having different numbers of dust effects active (sets C1-4), while Figure S12 compares cloud occurrence in these simulations. Based on cloud cover differences between simulation sets (see Figure S13), we estimate that the decomposition method could cause up to a 20% misattribution of forcing mechanisms in some regions. Notably, impacts of the partitioning method on low cloud cover ( Figure S13c) likely caused an underestimate of ERFs from warm cloud indirect effects over Africa (with compensating errors to ERFs attributed to other mechanisms). Our method was chosen among other possibilities to enable the sum of calculated dust-cloud effects to equal-with the exception of the small second "remaining" effect-the clean-air cloud forcing in the most standard simulation set (C1). Another option would have been to run all simulations with up to one effect removed at a time and compare, which would keep the simulations more realistic but potentially add to a different sum. (1)

Results
(d) Figure 2. Impacts of added mineral dust on clouds. Shown are simulated changes to cloud water path (liquid + ice water paths) (a), cloud cover (b), ice crystal number (c), and cloud droplet number (d) due to a hypothetical widespread increase in dustiness. Global mean changes are shown above each plot, as both absolute and percent changes from the present-day case. Clouds are separated between "low" (>700 mb), "mid-level" (400-700 mb), and "high" (<400 mb) levels, with (a,b) representing vertical integrations across each level and (c,d) representing averages among all ice and liquid clouds, respectively, in each level. Note that (c) represents all ice crystals rather than only those >5 μm as in Figure 1b. This figure was generated using output from simulation set C1.

Journal of Geophysical Research: Atmospheres
linkage, while Table 2 shows the associated global mean ERFs from all simulations. For simulations with our core estimates of input parameters (shown in Table 2 as bold), dust-cloud-climate interactions (semidirect plus all indirect) induced a + 0.02 ± 0.01 W/m 2 global mean ERF from an increase in dustiness and an ERF of −0.05 ± 0.02 W/m 2 for a decrease. We also performed sensitivity tests for each mechanism of dust-cloud-climate interactions. Combining the strongest cooling and warming responses across the four focused effects gives an ERF range from −0.26 to + 0.13 W/m 2 from a dust increase and

10.1029/2019JD031807
Journal of Geophysical Research: Atmospheres −0.21 to + 0.39 W/m 2 from a dust decrease (note that the most extreme forcings do not necessarily add linearly yet we present these summations because the simulations were designed to constrain only one dust-cloud effect at a time). Individual dust-cloud-climate links are described in this section, while direct and remaining effects-not the focus of this study-are shown in Figures S1 and S2 and briefly described in Texts S1 and S2, with Figure S3 supporting arguments therein.

Dust Indirect Effect via Cirrus Clouds
The addition of mineral dust in cirrus clouds had a mixed influence on ice crystal number density. While the additional source of INPs enables heterogeneous nucleation at relatively low supersaturations, this sink of water vapor prevents the air mass from reaching the high supersaturations necessary for homogeneous freezing of sulfate. Increasing dust led to an increase in ice crystal number over the North Africa region as well as Eurasia, which is downwind of both North African and Asian dust sources. Here, the new heterogeneous ice crystal source more than makes up for the lost source from homogeneous nucleation. Meanwhile, most other regions exhibit a decreased ice number due to the homogeneous effect dominating. This situation matches the spatial pattern of high cloud ice number changes in Figure 2c. Globally, the negative relationship between dustiness and cirrus ice crystal number dominates (Figure 1b). Note that increased INPs in mixed-phase clouds also contribute to the high cloud change, since a minority of mixed-phase clouds are above 400 mb. The cirrus ice crystal density impacts are associated with inverse effects on crystal size ( Figure S4) and proportional effects on cloud optical thickness.

Journal of Geophysical Research: Atmospheres
As shown in Table 2, the mineral dust increase (decrease) resulted in a net cirrus indirect effect ERF of −0.01 ± 0.02 (−0.05 ± 0.01) W/m 2 , with sensitivity tests demonstrating a possible range of −0.06 to + 0.01 (−0.05 to + 0.32) W/m 2 . As shown in Figure 4 for the case of increased dust levels, the concurrent effects on shortwave and longwave radiation have large global means dominated by changes in areas of cirrus thinning far from major deserts, but these predominantly cancel to render a weak global mean net effect. The net radiative effect roughly reflects the dust-induced changes in ice crystal density, with areas of increased (decreased) ice crystal number generally having a warming (cooling) net effect from the dust indirect effect via cirrus. Since cirrus clouds overall have a net warming effect on Earth's climate by trapping more longwave radiation than they reflect in the shortwave part of the spectrum, the small globally cooling net effect from increasing dust evidences a slight dominance of thinning cirrus due to the decline in homogeneous nucleation. Note that decreasing global dustiness has the opposite impacts on cirrus properties ( Figure 3c) and radiation ( Figure 5), but with larger differences per unit change in dust concentration. This is because the influence of dust on cirrus ice crystal density and radiative properties is strongest where homogeneous nucleation is the dominant ice crystal formation mechanism, as is generally the case with low global dust levels. Correspondingly, compared to the dust increase, the dust decrease produces a significantly smaller region where the dust change and high ice number impact are proportional (only over sections of Africa and Central Asia, as seen in Figure 3c). The dust decrease has the strongest impacts on longwave and shortwave radiation around 20°N while for the increase these are strongest around 10°S (see Figure S5). Our results could be biased by the model overestimating high cloud fraction in the equatorial and polar  Figure S14). The model also overestimates cloud radiative forcings in the tropics relative to Clouds and the Earth's Radiant Energy System -Energy Balanced And Filled (CERES-EBAF) v4.1 retrievals (Loeb et al., 2018) (see Figure S15). This overestimate is exacerbated in our model version relative to the bias of the default model (Bogenschutz et al., 2018). This may either be a result of our nudging to meteorology or a trade-off of the microphysical changes we implemented to more realistically represent dust-cloud interactions. As the net ERF results from a small remainder between complexly canceling shortwave and longwave components, it is difficult to estimate how results would be without these biases.
ERFs resulting from sensitivity tests are shown in Table 2. Lowering the deposition coefficient strengthened dust indirect effects via cirrus, while lowering or raising the updrafts fed into cirrus nucleation had a mix of results difficult to define globally. As expected, weakening dust INP efficiency as if it were aged uniformly reduced the associated forcings. As our model did not account for dynamic aging and hence we could only treat dust as if it were uniformly aged, our treatment here would only represent actual aging if this process quickly reaches completion. Cirrus indirect effects have opposite signs near and far from the major dust sources, and a more realistic aging process might weaken the latter with less impact on the former. Since indirect effects far from the major deserts dominate the global signal without aging, a more realistic aging could potentially weaken or even reverse the sign of the ERFs.
Use of the Dav19 PCF nucleation scheme in place of Ull17 resulted in substantially stronger longwave and shortwave impacts and affected their balance, such that a decrease in dust results in a moderately strong positive ERF of + 0.32 W/m 2 . Dav19 enables INP activation to occur where it would be negligibly weak in Ull17, producing cirrus with >10/L activated INPs from~100 mb lower in the atmosphere than with Ull17 (see Figure S7). Hence, using the Dav19 scheme enables altered dust concentrations to influence the properties and radiative effects of lower level cirrus that would not be significantly affected with Ull17. This is because the Dav19 scheme nucleates some ice at warmer temperatures and lower supersaturations than occurs in Ull17, as is evident when comparing the parameterization equations along with Note. Shortwave (SW), longwave (LW) and net (SW + LW) ERFs are shown. Simulations in bold represent core estimates while the rest represents sensitivity tests for factors that are not well constrained and simulations with alternate ice nucleation parameterizations. Included uncertainties are the standard deviation of global-mean impacts among the 5 years simulated, which did not represent interannual variability in dust emission.

Journal of Geophysical Research: Atmospheres
representative dust size values. In Ull17 the transition between weak and strong nucleation is sharper than in Dav19, which is likely an outcome of the INAS density formulation used in the Ull17 parameterization. The increased heterogeneous nucleation in Dav19 does not carry over to cirrus already having dense activated INPs with Ull17, and in cirrus with the highest supersaturations and lowest temperatures Dav19 is expected to produce less activated INPs than Ull17.

Dust Indirect Effect via Mixed-Phase Clouds
In mixed-phase clouds, an increase in INPs could be expected to primarily cause an increase in ice crystal number, since homogeneous nucleation does not compete in the mixed-phase temperature range like in cirrus. One caveat is that increased INPs could reduce ice crystal number in mixed-phase clouds by depleting cloud liquid earlier through the WBF process. This second effect has been found to be important when comparing a zero-INP atmosphere to one with dust INPs (Shi & Liu, 2019), though our simulations show that the positive relationship dominates during more natural INP perturbations. This is evident in how increasing the dust concentration has a positive impact on ice crystal number in midlevel clouds over a large region affected by North African and Asian dust (Figure 2c). The higher resulting ice crystal number is associated with smaller crystals (see Figure S4) and greater cloud optical thickness. Note that low and high clouds may also be affected by the dust change in mixed-phase clouds, as some of these clouds lie above 400 mb and below 900 mb (Shi & Liu, 2019). The positive relationship between dust and ice crystal number is visible in global mean values of ice number along mixed-phase cloud isotherms for the different dust levels ( Figure 1b). Dust changes are found to influence cloud phase, being negatively related to supercooled liquid fraction (Figure 1a). The cloud phase impacts found here are relatively modest, which may indicate that glaciation is largely enabled by ice from sources other than local ice nucleation (e.g., sedimentation of overlying cirrus crystals and ice detrainment from convective cores). While ice nucleation is generally increasingly efficient at colder temperatures, INPs and droplets become less available at higher altitudes. Hence, the influence of dust on ice number varies with temperature and has a maximum impact near −35°C (as seen in Figure 1b for glaciated clouds).
We find that enacting an increase (decrease) in dust concentrations causes a + 0.04 ± 0.01 (−0.07 ± 0.02) W/m 2 indirect effect ERF via mixed-phase clouds, with a + 0.01 to + 0.06 (−0.10 to −0.05) W/m 2 range from sensitivity tests. These clouds cause nearly equal shortwave cooling and longwave warming, and hence, changes to their optical thickness are expected to also have little net impact. Hence, while the shortwave and longwave ERFs due to changes in dustiness shown in Figures 4 and 5 and Table 2 are moderate, the effects are so balanced that the net effects (also shown) are small, with increased (decreased) dustiness globally having a slight warming (cooling) effect. Note that results may be biased due to a sizable overestimate of mixed-phase cloud cover. Figure S14 shows model comparisons to satellite retrievals of midlevel cloud cover, along with other levels that may contain mixed-phase clouds. Hence, longwave and shortwave mixed-phase indirect ERFs might be overestimated due to the cloud fraction bias. Model overestimates of tropical cloud radiative effects ( Figure S15) may also induce bias. It is uncertain how these factors could influence the calculated net ERFs, as biases to the occurrence and net radiative effects of affected mixed-phase clouds are difficult to distinguish from those of other clouds.
Inputting higher (lower) updrafts in mixed-phase clouds while increasing dust levels leads to increases (decreases) in ice nucleation rates and ice number density, but the longwave and shortwave effects still cancel to a slightly stronger warming (cooling) effect, with the opposite changes for a decrease in dustiness. Use of the DM15 immersion freezing parameterization in place of Ull17 induced weaker heterogeneous freezing rates overall (see Figure S11), leading to global effects of the same sign but weaker impacts on ice crystal density and radiation.

Dust Indirect Effect via Warm Clouds
Our simulations show increases in dust to slightly increase droplet number concentrations in low clouds around Africa (Figure 2d), which themselves have a generally cooling influence on Earth's radiative balance (shortwave scattering dominant). Dust affects droplet number in midlevel clouds as well, and feedbacks may exist with mixed-phase cloud indirect effects due to CCN effects influencing the number of droplets available for freezing. The associated warm cloud indirect effect ERFs shown in Figure 4 are dominated by a cooling effect strongest over the North Atlantic just west of the Sahel, an area with high occurrence of low clouds 10.1029/2019JD031807

Journal of Geophysical Research: Atmospheres
that become more reflective in reaction to the dust increase (the opposite case for a dust decrease is shown in Figure 5). However, dust CCN impacts are overall weak in our simulations. This is likely due to the atmosphere's abundance of CCN of other species, which limit dust's influence on droplet number by competing for a finite supply of water vapor. Dust is mostly insoluble and typically must mix with hygroscopic aerosol species to activate droplets. While aged dust particles of relatively large size can activate at low supersaturations, the ensuing condensation can prevent other CCN from forming droplets (Karydis et al., 2017).
Our calculations find a warm cloud indirect effect ERF from increased (decreased) dust of −0.01 ± 0.01 (−0.01 ± 0.01) W/m 2 , with sensitivity tests demonstrating a range of −0.02 to 0.00 (−0.01 to 0.00) W/m 2 . While this cooling due to low clouds appears dominant in the global net warm cloud indirect effect, effects in the longwave and shortwave each include a mix of cooling and warming effects between regions. As these signals tend to cancel between longwave and shortwave, they may indicate additional impacts of dust acting as CCN in higher liquid clouds, where cloud forcings are similarly balanced. Droplet number concentrations along with shortwave and longwave impacts were found to scale with the updraft speed fed into the aerosol activation scheme, with the net impacts of dust perturbations staying close to zero in all simulation sets. Though the calculated warm indirect effect ERFs are small, these may be biased toward higher magnitudes due to model cloud fractions being overestimated for extratropical low cloud, tropical midlevel clouds, and polar midlevel clouds ( Figure S14), along with tropical cloud forcings also being overestimated ( Figure S15).

Dust Semidirect Effect
The changes to low cloud cover shown in Figures 2b and 3b occur largely due to the altered temperatures and vertical temperature gradients resulting from dust absorption of shortwave radiation and consequent dust layer heating. Semidirect effects additionally influence cloud fraction in midlevel clouds (also shown in these figures), a region where mixed-phase indirect effects may also affect occurrence. Likewise, the longwave and shortwave semidirect responses to increased dustiness shown in Figure 4 show the influence of changes to both low and midlevel cloud cover. However, since the longwave and shortwave radiative effects of midlevel cloud cover responses largely cancel, the net semidirect effect primarily reflects differences in low cloud cover. In our simulations increased dust results in a greater low cloud cover over land areas downwind of the major dust sources, including much of Asia and Africa, while the reverse effect is evident over the Atlantic just west of the African continent. The decrease in low cloud cover narrowly outweighs the increase, and since low clouds have a net cooling effect the globally dominant radiative impact of increased dust is a net warming semidirect effect (Figure 4). For a large global decrease in dust levels, the opposite cloud and radiative changes occur.
We found a semidirect effect ERF from increased (decreased) dust of + 0.02 ± 0.02 (−0.02 ± 0.03) W/m 2 , while sensitivity tests give a range of −0.18 to + 0.06 (−0.04 to + 0.12) W/m 2 . The use of a higher dust imaginary refractive index for visible wavelengths, corresponding to a more absorbing dust layer, intensifies the ability for dust to raise low cloud fraction, switching the global mean net semidirect effect of a dust increase (decrease) to a cooling (warming). The calculated semidirect effects may be overestimated due to model biases in cloud forcing and fraction, for similar reasons as discussed above for warm cloud indirect effects.

Conclusions
In this study we estimated the differences in Earth's radiative balance (ERF) that could occur due to extreme changes in global dust levels through dust-cloud-climate interactions in the present-day atmosphere. Of the four mechanisms focused on (indirect via cirrus, mixed-phase, and warm clouds, plus semidirect effects), all showed weak global mean net radiative responses to dust emissions in the majority of simulation sets. All mechanisms responded to altered dust levels with a mix of warming and cooling effects among regions for both the longwave and shortwave TOA energy balances, which themselves largely canceled to render smaller net impacts.
We found that increasing dust is most likely to cause a slight global cooling through cirrus indirect effects. This conclusion is in line with the findings of Gasparini and Lohmann (2016) that global impacts of added INPs are weak due to compensating influences on ice crystal number and size. By contrast Storelvmo et al. (2013) found that the addition of INPs to cirrus clouds could cause substantial global cooling. However, in that study INP concentrations were held within a specific range globally, as is unlikely to occur for a natural aerosol like dust that varies in abundance over both time and distance from its emission source. Other studies have found strong cooling effects due to ice nucleation on dust in cirrus (Kuebbeler et al., 2014;Liu et al., 2012), but these compared cases with and without dust INPs rather than dust perturbations from realistic levels.
Increased dust caused a slight warming through indirect effects in mixed-phase clouds in our simulations. We also found dust INP abundance to be associated with optically thicker mixed-phase clouds. This contrasts with the findings of Shi and Liu (2019), though that study assessed the impact of all dust INPs collectively, for which the ability of INPs to induce cloud liquid depletion was likely more prominent than for perturbations between more natural states. The ERFs reported here are both weaker and of the opposite sign as previously concluded in Sagoo and Storelvmo (2017). This is likely related to our simulations having a different implementation of immersion freezing that resulted in less ice nucleation in clouds warm enough to have a net cooling influence. Additionally, Sagoo and Storelvmo (2017) simulated a pre-industrial climate rather than the present-day atmosphere considered here. Anthropogenic aerosol acting as INPs may render indirect effects weaker than in past climates, though our simulations did not include these species as INPs.
In warm clouds, increased dustiness was found to slightly enhance cloud droplet concentrations, with an ensuing small cooling around and near North African deserts, and the opposite case for decreased dustiness. Our simulations show the opposite semidirect effects as observations over the tropical North Atlantic just west of Africa (Amiri-Farahani et al., 2017;Doherty & Evan, 2014) and a modeling study (DeFlorio et al., 2014). This may be linked to our use of fixed SSTs, which prohibited the ocean surface from cooling in response to overhead dust blocking sunlight. The present-day atmosphere contains abundant anthropogenic aerosols that limit dust's influence on warm clouds by acting as a stable supply of CCN and masking over some of the light otherwise scattered and absorbed by dust.
While the majority of simulation sets herein found weak global net radiative impacts of changes to dust emissions through clouds, moderate regional sensitivity was found in areas with clouds near the major dust sources. Considering how greatly the studied net ERFs varied from region-to-region, dust changes to specific desert areas could potentially have a greater global impact than the globally uniform changes simulated here. Also, the global shortwave and longwave radiation budgets were found to be sensitive to dustiness, suggesting that dust may have a stronger effect in global climate simulations that allow the atmosphere and oceans to fully adjust to the dust perturbations.

Data Availability Statement
Data from the CESM simulations described herein are shown in Table 2 and presented in figures both in the text and supporting information. Output files from each of these simulations and dust emission files used are available online (https://doi.org/10.5281/zenodo.3701671) for selected variables as 5-year averages. Model setup and source code modifications required to replicate these simulations are described in sections 2.1-2.3, while analysis methods are detailed in section 2.4.