Potential limitations of using a modal aerosol approach for sulfate geoengineering applications in climate models

geoengineering applications in climate models Daniele Visioni1, Simone Tilmes2, Charles Bardeen2, Michael Mills2, Douglas G. MacMartin1, Ben Kravitz3,4, and Jadwiga Richter5 1Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, NY, USA 2Atmospheric Chemistry, Observations, and Modeling Laboratory, National Center for Atmospheric Research, Boulder CO, USA 3Department of Earth and Atmospheric Sciences, Indiana University, Bloomington, IN, USA 4Atmospheric Sciences and Global Change Division, Pacific Northwest National Laboratory, Richland, WA, USA 5Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder CO, USA Correspondence: Daniele Visioni (daniele.visioni@cornell.edu)

probability distribution from Kärcher and Burkhardt (2008) that determines when the supersaturation is above the threshold for homogeneous freezing to happen. This means that homogeneous freezing can happen only when local vertical velocities are high (Kärcher and Lohmann (2002)). Those velocities are determined following Morrison and Pinto (2005) as with the turbulent kinetic energy (TKE) determined using a steady state energy balance as in Bretherton and Park (2009).
For heterogeneous freezing, only Coarse mode dust particles are assumed to be available IN (although the exclusion of soot particles and other particles as inefficient ice nuclei has been debated; see Kärcher and Burkhardt (2008)). Given the internal mixing assumption in MAM3, the number of available dust nuclei is determined as a fraction of the overall amount of Coarse 165 mode particles given the mass of dust and the overall aerosol mass where m d is the mass of dust in the Coarse mode, m ss is the mass of sea salt in the Coarse mode, and N c is the total number of particles in the Coarse mode. This approach assumes a negligible amount of Coarse mode sulfate in the upper tropical stratosphere in the denominator of the fraction. Another shortcoming of the heterogeneous freezing description in this version 170 of WACCM, as already discussed by Mills et al. (2017), is that when heterogeneous freezing occurs, IN that have nucleated to form ice particles are not removed from the available reservoir, thus allowing for too many particles to be formed via heterogeneous freezing.
This approach to the microphysical modeling of cirrus clouds in CESM has been discussed recently by Maloney et al. (2019), comparing it with a more complex approach using the coupled Community Aerosol and Radiation Model for Atmospheres by Junge et al. (1961). This happens because of the presence of various sources of sulfate aerosols in the stratosphere, even without considering volcanic aerosols from explosive volcanic eruptions. In particular, surface emissions of carbonyl sulfide (OCS) and dimethyl sulfide (DMS), which are light, well mixed gases that may reach the stratosphere where they are oxidized, forming sulfate aerosols (Vet et al. (2014)). Meteoric sulfur also plays a part in the formation of the layer (see Gómez Martín et al. (2017)). The quantity of sulfate aerosols produced in the stratosphere with an artificial injection of SO 2 would be larger 195 by some orders of magnitude compared to the amount in the quiescent Junge layer; the full stratospheric distribution in the case of GLENS has been described elsewhere (i.e., Kravitz et al. (2019)), so here we simply report that at 50 hPa, between 30 • N and 30 • S, the average mass concentration for the last 20 years of simulation is 163 µg/kg, compared to 0.4 µg/kg in the unperturbed case; while this simulation considered high cooling (∼ 4 • C reduction in global mean temperature), even much smaller amounts of SAI would lead to substantially more stratospheric sulfate than the unperturbed case. This magnitude is 200 driven by the long lifetime of the produced aerosol particles in the stratosphere (around 12 months, Visioni et al. (2018b)).
Once the particles cross the tropopause, the various removal processes strongly reduce the lifetime, driving the concentration down. This is visible already in Fig. 1a, with the tropopause clearly visible as a change in concentration. Nonetheless, the average concentration of sulfate under SAI is much larger even in the upper troposphere ( Fig. 1b) and only returns close to the Baseline levels close to the surface, where other sulfate sources are predominant (Fig. 1c).
205 Fig. 1 shows the behavior of dust and sea salt in the Coarse mode in the same altitude-latitude region (those are the only other aerosol species considered in this version of CESM1 in the Coarse mode: black carbon is only found in the Accumulation mode). Various behaviors can be observed depending on the altitude of analyses for the GLENS case. We focus here on the uppermost troposphere, right below the tropopause (black boxes in Fig. 1) and on on a layer below, in the upper troposphere (green boxes). In the upper level, both dust and sea salt present a behavior that is much more similar to that of sulfate: a large, 210 initial increase in concentration in the first 5 to 10 years followed by either a constant evolution or by a small decrease (from 0.1 to 0.2 µm/kg for dust in the upper layer, from 0.06 to 0.2 µm/kg and from 0.9 to 0.17 µm/kg for sea salt in the upper and lower layer). Only dust in the lower layer does not show such an abrupt initial increase. On the other hand, in the Baseline case the dust mixing ratios show a decrease by 25% over the 80 years of analyses in the lower level. There is no immediate physical mechanism by which such changes would be observed in aerosol species independent of sulfate; therefore, a more in-depth 215 analysis is necessary.
One possible explanation would be a change in the surface emissions of these species. However, the vertical profiles of tropical concentration for the three species seem to exclude that. This is further confirmed by the analyses of the overall burdens shown in Fig. 2. For these, no initial abrupt change can be observed, and the time evolution shows the opposite behavior from that observed in Fig. 1, with a minor reduction in GLENS compared to Baseline. The behavior of the Baseline case is 220 consistent with previous projections: Mahowald and Luo (2003) predicted a reduction in overall dust emission with increasing GHGs due to higher precipitation and more surface moisture produced by the warmer air; thus the precipitation decrease in GLENS (Cheng et al. (2019)) would exacerbate differences between dust burden in GLENS and the baseline. Struthers et al. (2013) on the other hand projected a small increase in sea salt aerosols mostly as a result of increased surface wind speeds.
Having excluded that changes in surface emissions are the cause of our observed abrupt increase, we further investigate the 225 smaller aerosol modes (Fig. 3). In this case, we see a much closer agreement of the values in GLENS to those of the initial baseline conditions for sea salt and dust: this indicates that the initial aim of the GLENS simulations, which was to maintain the state of the climate as close as possible to 2010-2029 conditions, also helps with not modifying surface emission sources of those species. The time evolution of the globally averaged quantities for these species and modes also doesn't show the same abrupt change as that observed in Fig. 1 (see Fig. S1 and S2). So clearly the cause of the rapid change in concentration of 230 the Coarse mode must be the result of the very rapid increase in SO 4 aerosol concentration by more than 4 orders of magnitude.
The solution to this conundrum can be found by analyzing the behavior of the simulated radius of the particles in the three modes ( Fig. 4). In MAM3, all aerosol species are assumed to be internally mixed. The simplified assumption that different aerosols peak at different locations in the lower atmosphere is reasonable for the baseline scenario ( Fig. 1 for the vertical 235 distribution, and Fig. 2 for the spatial distribution). The parameterization requires that in each gridbox, all aerosols in one mode are treated as the same in terms of their size distribution, so for each mode, only one mode diameter and one number concentration is used in the online calculations (since the geometric standard deviation is constant throughout the atmosphere).
The mass concentrations of the different species are then calculated using a reference density for each species (Liu et al. (2012)). This information is necessary to calculate the gravitational settling velocities (at each layer), which follow Seinfeld 240 and Pandis (2016), where the equation of a free-falling spherical particle in a fluid which has reached terminal velocity (which is done in less than 1 × 10 −6 for particles of diameter of 1 µm) is shown to be: where D p is the diameter of the particle, ρ p is the density of the particle, C c is the slip-correction factor that accounts for non-continuum effects (dependent on the diameter of the particles), and C D is the empirical drag coefficient (dependent on the 245 Reynolds number R e ).
The initial observed changes can thus be explained this way. The SO 4 formed in the stratosphere is the predominant form of aerosols in the stratosphere, deciding the radius of the particles at those altitudes. Most of these particles are in the Coarse mode, therefore changing the radius considered by the model in those grid-boxes. Once the sulfate aerosols cross the tropopause (whether by gravitational settling, mostly in the tropics, or by large scale circulation, at higher latitudes) they are added to the 250 particle distribution already present in the upper troposphere. Due to the small amount of Coarse mode aerosols there, and due to the internal mixing assumption in MAM3, the dust and sea salt aerosols already present are "forced" to reduce in size (and to conserve mass, to increase in number concentration, see Fig. S3). This produces a drop in the overall size by around 25% within 3 years in the upper troposphere. Because of the strong dependence of gravitational settling velocities on size, the result is a drop in these velocities that, for dust aerosols, can be estimated to be over 40%. This is not a direct output of the model,    by decreasing removal, especially in regions where contributions from below are scarce. And the reason why this is only visible in the upper troposphere is that, further below, the preponderance of pre-existing, background Coarse mode aerosol reduces the strength of this effect, not changing the radii. The preponderance of other removal mechanisms (in-cloud scavenging) also We have shown that the simulated changes in dust and sea salt in the upper troposphere in GLENS using CESM1(WACCM), are non-physical in nature. However, by looking at the concentration of the same species in the Baseline conditions it is clear that the overall presence of those aerosols at the altitudes of analyses is small, and therefore would have a negligible effect both 265 when looking at the direct effects (for instance, on the radiative fluxes, which would be overshadowed by the effects of the sulfate aerosols in the stratosphere) and when eventually looking at surface effects. We show in this section, however, that their change does influence the simulation of ice clouds formation in the upper troposphere. Visioni et al. (2021) showed that some noticeable changes in cloud cover are present in the GLENS simulations, and that when separating the contribution of different types of clouds, most of these changes were attributable to changes in high clouds, that are defined in CESM1(WACCM) as all 270 clouds formed between the altitude of 400 and 50 hPa. Furthermore, these changes were not as noticeable in other performed simulations where the same cooling as GLENS was achieved using a solar constant reduction approach, and when together with that the stratospheric heating produced by the aerosols was imposed, without the presence of the aerosol themselves. This would point towards a contribution to the observed changes in high cloud cover in GLENS by some effects connected to the aerosol themselves. In particular, this points to some changes in the freezing processes that produce ice crystals in the upper 275 troposphere, given that at those altitudes there is much less water than in the lower troposphere, and most of that is in the form of ice rather than liquid droplets (see Fig. S5).
As we detailed in Section 2.2, there are two types of processes that may lead to the formation of ice crystals at those altitudes: the spontaneous freezing of small aqueous sulfate droplets at high relative humidity rates (Chen et al. (2000)), known 280 as homogeneous freezing, and the freezing of water droplets mediated by the presence of insoluble aerosol particles (ice nuclei, IN), known as immersion or heterogeneous freezing (Diehl and Wurzler (2004)), which can happen at lower relative humidity conditions and higher temperatures (Knopf and Koop (2006)). Due to the difficulties in measuring the amount of ice crystals in the upper troposphere, and to the challenges of representing the processes in models, there are plenty of uncertainties over the predominance of one of the two formation processes over the other. In CESM, only sulfate particles in the Aitken mode can 285 act as the substrate over which homogeneous freezing can take place; this is generally understood to be correct, as too large particles would be increasingly harder to freeze (Chen et al. (2000)). Since Aitken mode particles do not change in GLENS (see the model "sees" is therefore not physical in origin, but simply an artifact of the microphysical parameterization and thus should not be treated as a physical side-effect to SAI. A clue as to the two different mechanisms at play can be determined by observing the rate of change at low versus high latitudes. At low latitudes, the changes happen gradually as the stratospheric sulfate load (and thus the warming) increases. At high altitudes, the in-cloud ice changes are very abrupt, much more similar to the aerosol changes observed previously in this work. Further confirmation can be derived by observing the differences 310 between the full GLENS simulations and the simulations described by Visioni et al. (2021) with stratospheric heating imposed, but no aerosols. By comparing the two (see bottom panels of Fig. 5), the effect of the presence of the aerosols can be separated from the dynamical ones, and indeed this confirms the different nature of the observed ice changes, since the two simulations (aerosols+stratospheric heating and stratospheric heating alone) show no differences at low latitudes, whereas they are as different as GLENS-Baseline at high latitudes.

315
Further information would be gained by observing the changes in heterogeneous versus homogeneous freezing in these simulations. However, output fields that separate the two processes have not been saved in the original simulations from GLENS. Therefore, a complete analysis is not possible using the whole ensemble and time period. Noting that the main, unexplained processes (with our current understanding of ice formation) happen in the very first decade of the simulations, for the following analyses we have re-run the first 21 years of one of the ensemble members with the same configuration (thus 320 bit-by-bit, the result is the same) but retaining information related to the two different upper tropospheric freezing processes.
The results in Fig. 6 confirm our previous observations. For homogeneous freezing, very few changes are initially observable at all latitudes, whereas large changes are observable in the concentration of what the model considers to be IN for heterogeneous freezing, sharply in the first few years. This increase is much larger than the increase in Coarse mode dust aerosols observed in Fig. 1, where the overall concentration was doubled, while here the increase is much larger, at more than 10 times the amount 325 in the first two years. To partially explain this change we refer to the equation used by the model to determine the fraction of the available dust in the Coarse mode, which is used as IN for heterogeneous freezing (Eq. 2.2). The available number of dust particles that can be used by the model depends on the overall amount of Coarse mode particles, but is weighted by which fraction, locally, of aerosol is dust compared to sea salt. This formula (which in normal conditions gives satisfactorily results, see Maloney et al. (2019)) presumes a lack of sulfate in the upper troposphere. In our case however, as we have shown, the amount of Coarse mode sulfate is actually by far preponderant in the upper troposphere. So while the number of Coarse mode particles N C grows (see Fig. S3), the fraction of all aerosols that is considered does not. This is a separated, but similar, problem as the one observed for the behavior of the aerosols. We tested how sensitive this assumption is by performing an identical simulation to the one that gave us information about homogeneous and heterogeneous freezing, but modifying the code so that Eq. 2.2 becomes in the first 21 years of simulation, for the three latitudinal band already considered in Fig. 5. On the right, the column mean between 50 and 400 hPa is considered (thick lines), together with the case described in the text with a sensitivity test for Eq. 2.2 (dotted lines).In panel f) the ice number concentration shown in Fig. 5 is compared with that for the sensitivity experiment.
Here we investigate whether the changes shown in the previous section produce some signal in the modeled radiative fluxes at the top of the atmosphere, as that would be important in order to determine how significant the effect is. To do so we use the method suggested by Ghan (2013) and the forcing produced by the changes in clouds (∆C) can be defined as As Ghan (2013) notes, normally these fluxes are just an estimation of the real components of the forcing, as to properly estimate 355 them it would be necessary to maintain surface and tropospheric temperatures fixed between the two cases to avoid changes in the radiative emissions of the troposphere, which would be at different temperature in the two cases. In our case, however, the GLENS simulations have been performed on purpose to maintain similar global surface temperatures as the Baseline 2010-2030 period: some small regional temperature differences are still present ), but overall this can be considered a minor factor and the estimated forcing quite robust if we are comparing the radiative fluxes always against the 360 2010-2030 period. In the case of cloud changes, this further assures that we are not counting effects produced by the surface warming, which might also locally modify cloud cover.
In Fig. 7 we show the global evolution of ∆F D and ∆C in GLENS. The aerosol direct radiative effect is linear, positive in the shortwave (implying a cooling) as they reflect incoming sunlight and negative in the longwave (implying a warming) as they absorb and re-emit longwave radiation, with an overall negative effect, as expected. We show the latitudinal breakdown 365 of the fluxes in panels 7b-c normalized by the stratospheric AOD produced by the SO 2 injections for the initial period (we pick 2026-2035 to avoid the very first few years when the algorithm that determines the SO 2 injections is still converging, see Kravitz et al. (2017)) and the last twenty years (2081-2100). A partial non-linearity can be observed between the two periods. This can be explained by the slight increase in stratospheric sulfate aerosol during the whole period of analyses (see Visioni et al. (2020a), the effective radius in the stratosphere in GLENS grows from 0.4 µm to over 0.5 µm as a result of the 370 increasing injection rates resulting in more coagulation of SO 2 with pre-existing particles). Larger aerosols scatter slightly less efficiently, but they absorb more LW radiation. The changes in upper tropospheric aerosols described before do not influence these radiative fluxes, as they are, in proportion, negligible compared to the stratospheric increase. Looking at the ∆C fluxes, however, we can see that the behavior of the LW ∆C shows a very different behavior compared with the SW. This effect is positive in the first 30 years, and then negative afterwards, whereas the SW forcing is always positive. The SW forcing is 375 easily explainable as directly connected to the presence of the aerosols above: if less sunlight reaches the troposphere, because it is partially reflected in the stratosphere, the same clouds would appear as to be less reflective, hence the positive sign of the SW (cloud masking). This doesn't apply to the LW, as globally the surface temperatures are the same, but can be explained by noting that the main contributor to LW trapping by clouds is upper tropospheric ice (Fusina et al. (2007)). Since, in the first decades, mid-latitudinal ice clouds increase as shown in Fig. 5 because of "more" dust aerosols acting as IN, they would 380 trap more outgoing LW radiation. When tropical ice clouds start decreasing because of the dynamical mechanism produced by stratospheric heating, the positive bump is erased by the negative forcing (less ice clouds, more LW radiation escaping to space) already analyzed by Visioni et al. (2018a). This is further confirmed by the latitudinal breakdown of ∆C in panels 7e-f, where the (normalized by stratospheric AOD) forcing in the LW goes from positive, in the first 10 years, to negative in the last 20.

Conclusions
In this work we have identified the presence of some weaknesses in the 3-mode modal approach (MAM3) used in CESM1(WACCM) when a large presence of aerosols settles down from the stratosphere, as under stratospheric sulfate aerosol injection, which results in some artificial changes to cirrus clouds. MAM3 only separates the species in three modes depending on their size by 390 treating all aerosol species as the same. When a large amount of sulfate is produced (mostly in the coarse, largest mode) in the stratosphere following the injection of SO 2 in large quantities, the particles slowly descend upon the troposphere, where they are quickly removed. However, the size of the Coarse mode particles is different (smaller) compared to that of the Coarse mode particles already present in the troposphere, whose source is at the surface: therefore, MAM3 "sees" an abrupt decrease in all aerosol size in each gridbox in the upper troposphere for the coarse mode, resulting in smaller settling velocities for aerosol 395 species that would not be otherwise affected in the real world. This effect results in an increase in the mass of particles of all species at those altitudes even if there would be no natural causes for it to happen. The effect is small, and its direct effect on the cooling produced by the stratospheric aerosols negligible. However, the unnatural addition of dust in the upper troposphere results in more particles that the model can use as solid ice nuclei for the freezing of ice particles in clouds. This effect is particularly evident at mid-and-high latitudes, where the low relative humidity and lack of aerosols make other ice formation 400 processes more scarce in normal conditions. The formation of these ice crystals, as simulated by the model, does produce a noticeable change in the cloud forcing, especially in the first years, when the effect from the incorrect presence of the dust aerosols is large compared to other, dynamical effects that tend to reduce the amount of ice crystals in the upper troposphere when a warming of the stratosphere is present. Given the set-up of the GLENS experiments, this effect is counteracted by the presence of the algorithm which determines how much SO 2 is needed every year to counteract the effect of the increased emissions in RCP8.5: even if the radiative imbalance were to be large, the algorithm would just prescribe more SO 2 to be injected, therefore resulting, overall, in the same global mean temperatures as if the error in the ice cloud formation was not there. In the sensitivity simulations produced for Section 4, which reduced the amount of ice clouds incorrectly formed by the model, the cumulative amount of SO 2 injected 410 in the first 20 years was 131 Tg-SO 2 (6.2 Tg-SO 2 /year), whereas the same period in the default simulation had a cumulative amount of SO 2 of 154 Tg-SO 2 (7.3 Tg-SO 2 /year): considering that, aside from the change in Eq. 2.2, the two simulations were otherwise the same, we can assume that this difference arises from the changes in ice clouds as simulated by the model.
Overall then, if we assume this effect is constant throughout the whole simulation period, it would account for an cumulative injected 88 Tg-SO 2 , in the entire century: while this amount might seem large, it accounts for less than two years of anthro-415 pogenic emissions of SO 2 at present (Visioni et al. (2020b)). This, furthermore, assumes a very large injection to overcome a considerable amount of warming (over 4 • ) over the entire century. More moderate mitigation scenarios would require far less cumulative injections of SO 2 (for instance, see Tilmes et al. (2020)).
This doesn't imply that the effect can be ignored, but suggests that going forward when using a modal approach to aerosol 420 microphysics, simulations where large amount of sulfate is present in the stratosphere should treat sulfate aerosols in a separate coarse mode that is not internally mixed with the other species compared (dust, sea salt, black carbon, etc.). There could be various applications where this observation might be useful: SAI is one example, but the simulation of explosive volcanic eruptions is another case where it could be useful: for instance, Schmidt et al. (2018) used the same CESM1(WACCM) model described here to estimate the global volcanic radiative forcing in the last 45 years, and ran into a similar observation as found in this 425 paper, but did not find an explanation for it. The mechanism that produced, in their simulations, an increase in ice particles in the upper troposphere is definitely the same we have encountered here, and would explain some of the LW forcing changes diagnosed in their simulations with a similar method (see Fig. 4 and 5, Schmidt et al. (2018)). More generally, it would be crucial to properly represent the upper troposphere in the case of volcanic eruptions to verify the influence of volcanic eruptions on ice clouds as observed by some studies in the observational record (see for instance Friberg et al. (2015)). Future modeling efforts models used for simulating sulfate geoengineering application are not applied outside of the parameters for which the models will give reliable answers. Most of the time, models are also compared against volcanic events to make sure they are properly simulating the results, but not all effects are immediately visible in those cases as compared to a long term geoengineering simulation, and a comparison of cloud changes after a volcanic eruption may be complicated by various factors (Friberg et al. (2015)).