Climate response to off-equatorial stratospheric sulfur injections in three Earth system models – Part 1: Experimental protocols and surface changes

There is now substantial literature on climate model studies of equatorial or tropical stratospheric SO₂ injections that aim to counteract the surface warming produced by rising concentrations of greenhouse gases. Here we present the results from the first systematic intercomparison of climate responses in three Earth system models wherein the injection of SO₂ occurs at different latitudes in the lower stratosphere: CESM2-WACCM6, UKESM1.0 and GISS-E2.1-G. The first two use a modal aerosol microphysics scheme, while two versions of GISS-E2.1-G use a bulk aerosol (One-Moment Aerosol, OMA) and a two-moment (Multiconfiguration Aerosol TRacker of mIXing state, MATRIX) microphysics approach, respectively. Our aim in this work is to determine commonalities and differences between the climate model responses in terms of the distribution of the optically reflective sulfate aerosols produced from the oxidation of SO₂ and in terms of the surface response to the resulting reduction in solar radiation. A focus on understanding the contribution of characteristics of models transport alongside their microphysical and chemical schemes, and on evaluating the resulting stratospheric responses in different models, is given in the companion paper (Bednarz et al., 2023). The goal of this exercise is not to evaluate these single-point injection simulations as stand-alone proposed strategies to counteract global warming; instead we determine sources and areas of agreement and uncertainty in the simulated responses and, ultimately, the possibility of designing a comprehensive intervention strategy capable of managing multiple simultaneous climate goals through the combination of different injection locations.

We find large disagreements between GISS-E2.1-G and the CESM2-WACCM6 and UKESM1.0 models regarding the magnitude of cooling per unit of aerosol optical depth (AOD) produced, which varies from 4.7 K per unit of AOD in CESM2-WACCM6 to 16.7 K in the GISS-E2.1-G version with two-moment aerosol microphysics. By normalizing the results with the global mean response in each of the models and thus assuming that the amount of SO₂ injected is a free parameter that can be managed independently, we highlight some commonalities in the overall distributions of the aerosols, in the inter-hemispheric surface temperature response and in shifts to the Intertropical Convergence Zone, as well as some areas of disagreement, such as the extent of the aerosol confinement in the equatorial region and the efficiency of the transport to polar latitudes. In conclusion, we demonstrate that it is possible to use these simulations to produce more comprehensive injection strategies in multiple climate models. However, large differences in the injection magnitudes can be expected, potentially increasing inter-model spreads in some stratospheric quantities (such as aerosol distribution) while reducing the spread in the surface response in terms of temperature and precipitation; furthermore, the selection of the injection locations may be dependent on the models' specific stratospheric transport.

1 Introduction

A climate model is an imperfect representation of the real climate system, having to deal with simplified processes across spatial and temporal scales, bound by both technical limits to the computational capabilities and fundamental uncertainties in some of the underlying processes. Reducing such uncertainties can be a complicated process and may leverage various strategies: (i) by simulating current and past climates for which we have available observations, the skills of current climate models can be evaluated and the reliability of future projections assessed (Brunner et al., 2020). (ii) The use of large ensembles of simulations can further help to narrow down some of the uncertainties by isolating biases due to internal model variability over time (Deser et al., 2012; Lehner et al., 2020). (iii) Finally, the use of multiple climate models can help to clarify some of the sources of uncertainty by analyzing the response to an external forcing in multiple independent models. This has been done in numerous generations of the Climate Model Intercomparison Project (CMIP) (Masson and Knutti, 2011; Boucher et al., 2013; Tebaldi et al., 2021), which addresses the effects of an increased CO₂ forcing on the global and regional surface climate (Hawkins and Sutton, 2011). Other modeling exercises such as the Chemistry–Climate Model Initiative (Morgenstern et al., 2018) have used models with explicit atmospheric chemistry coupling (which many of the CMIP models lack) to determine the future evolution of key atmospheric variables, such as stratospheric ozone (Dhomse et al., 2018) and stratospheric circulation (Eichinger et al., 2019). Finally, model intercomparisons have been designed to constrain the response of climate models to explosive volcanic eruptions, trying to understand both the evolution of the aerosol cloud over time (Timmreck et al., 2018) and its climatic effects (Zanchettin et al., 2016; Clyne et al., 2021; Zanchettin et al., 2022).

The assessment of simulated climate intervention techniques presents all of the challenges discussed above, but also novel ones. The proposed injection of sulfur dioxide (SO₂) to temporarily cool the planet (stratospheric aerosol injection, SAI) while greenhouse gas levels are reduced (Crutzen, 2006) has been studied with climate models in the last 2 decades, and multiple intercomparisons have been carried out spearheaded by the Geoengineering Model Intercomparison Project (GeoMIP, Kravitz et al., 2011). The first SAI experiments, GeoMIP G3 and G4 (Kravitz et al., 2011), prescribed the injection of SO₂ in the equatorial stratosphere (in a time-varying or constant manner, respectively). However, the modeling teams were free to specify the exact details of the injection “in the same manner as they simulate a volcanic eruption”; this led to differences in the approaches chosen. For instance, the University of L'Aquila Chemistry–Climate Model (ULAQ-CCM) injected SO₂ between 18 and 25 km of altitude using a Gaussian distribution, whereas the Goddard Chemistry–Climate Model (GEOSCCM) injected SO₂ uniformly between 16 and 25 km of altitude. The varying experimental protocols resulted in large differences in the simulated aerosol distributions, which compounded the differences coming from the differences in the representation of aerosol microphysics amongst the models (Visioni et al., 2017). More recently, the GeoMIP experiment G6 (Kravitz et al., 2015) prescribed variable injections of SO₂ between 10^∘ N and 10^∘ S and between 18 and 20 km of altitude. However, some models lacking interactive aerosol treatment prescribed an aerosol distribution from sources not entirely consistent with the original prescription for models without interactive microphysics, resulting in challenges in the attribution of the causes of some of the resulting discrepancies (Visioni et al., 2021). Lastly, the test-bed experiment described by Weisenstein et al. (2022) prescribed two injection scenarios (one uniformly between 18 and 20 km and between 30^∘ N and 30^∘ S; one in two precise locations, 30^∘ N and 30^∘ S, at 19 km of altitude) to understand inter-model differences between injections of SO₂ and direct injections of accumulation-mode H₂SO₄ to better control the size of the resulting aerosols.

In parallel, simulations with the Community Earth System Model version 1 with the Whole Atmosphere Community Climate Model (CESM1(WACCM)) were carried out to examine how changes in the location of the sulfate injection matter for the climatic impacts. Tilmes et al. (2017) injected SO₂ at four locations (30, 15^∘ N, 15 and 30^∘ S) and showed that they resulted in different aerosol distributions and in different atmospheric and surface responses (see also Richter et al., 2017; Tilmes et al., 2018b). Further, MacMartin et al. (2017) showed that these SO₂ injection locations could be combined to achieve multiple climate goals with the use of a feedback algorithm capable of deciding each year how much SO₂ to inject at each of the aforementioned latitudes (Kravitz et al., 2017). The four locations of injection were necessary to control the overall geographical distribution of the simulated aerosols and thus to constrain not only the global mean surface temperatures, but also to make sure that the large-scale inter-hemispheric surface temperature gradient and the Equator-to-pole surface temperature gradient are maintained. This approach was adopted in the Geoengineering Large Ensemble experiment (GLENS, Tilmes et al., 2018a) that studied the feasibility and climate impacts of a potential SAI deployment strategy. By differentiating the injection locations, numerous side effects and adverse impacts of SAI could be reduced compared to an equatorial injection strategy (Kravitz et al., 2019). Up to now, all of the above assessments and the exploration of other strategies have been performed with only one climate model – CESM1(WACCM).

In this work we aim to perform the first systematic intercomparison of the stratospheric and surface climate responses to the injection of SO₂ at different latitudes in the stratosphere. We use three comprehensive Earth system models (ESMs) with interactive sulfate aerosol treatment and a set of carefully designed experiments with a single point of injection per experiment. Our goals are (i) to robustly evaluate similarities and differences in the simulated sulfate aerosol distributions and in the resulting surface and atmospheric responses; (ii) to elucidate the areas and sources of inter-model differences, in particular the roles of the specific model microphysical schemes and the biases in climatological transport; (iii) to determine the reliability of estimates of the surface responses to SAI previously performed in MacMartin et al. (2017) as well as the middle-atmospheric responses performed in Tilmes et al. (2018b) and Richter et al. (2017) using the multi-model approach; and (iv) to lay the groundwork for an inter-model comparison of SAI simulations achieving the same set of surface temperature goals defined in Kravitz et al. (2017) and mentioned above using a feedback algorithm. In this paper we discuss the simulated aerosol fields and their effects on zonal mean surface temperatures and precipitation; we then use that information to determine the magnitude of SO₂ injection needed in each location to obtain the desired temperature targets. In the companion paper (Bednarz et al., 2023, hereafter Part 2), the simulated differences in the aerosol distribution are explained in terms of the model differences in atmospheric circulation, and the resulting SAI impacts on stratospheric temperatures, chemistry and circulation are discussed.

2 Earth system models used

We performed our simulations with three Earth system models (ESMs): CESM2-WACCM6, UKESM1.0 and GISS-E2.1-G; for the latter, two different versions with two different aerosol treatments were used. One simulation for each model has been performed. All models have full atmosphere–ocean–land coupling as well as an explicit aerosol treatment and stratospheric chemistry. A brief description is provided below for all of the models, and a discussion of the differences in the aerosol treatments of particular importance for this work is given at the end of this section.

2.1 CESM2-WACCM6

The Community Earth System Model version 2 with the Whole Atmosphere Community Climate Model version 6 (CESM2-WACCM6, Gettelman et al., 2019; Danabasoglu et al., 2020, hereafter CESM2) is used with comprehensive stratospheric and upper-atmospheric chemistry as well as interactive aerosol microphysics using the Modal Aerosol Module (MAM4) (Liu et al., 2016). Its horizontal resolution is 1.25^∘ longitude by 0.9^∘ latitude, with 70 vertical levels and a model top at about 140 km. An evaluation of the model response to past volcanic eruptions was performed in Mills et al. (2017) and Schmidt et al. (2018) using the previous version of the model, CESM1(WACCM). Unlike the CESM2 version described in Danabasoglu et al. (2020), which includes both comprehensive tropospheric and stratospheric chemistry, the version used here includes comprehensive stratospheric chemistry but only a simplified chemistry of importance in the troposphere; this is thus similar to the chemistry scheme in the CESM1(WACCM) version used for previous geoengineering studies.

2.2 GISS-E2.1-G

NASA's ModelE is developed by the Goddard Institute for Space Studies (GISS). The current version, ModelE2.1, has a horizontal resolution of 2.5^∘ longitude by 2^∘ latitude and 40 vertical layers extending up through the mesosphere (model top of ∼80 km). The version used in this study, GISS-E2.1-G, has a fully interactive ocean with 32 vertical layers based on the Russell ocean (Kelley et al., 2020). While GISS ModelE also has the ability to use the HYCOM ocean (Kelley et al., 2020), that option is not employed here. GISS has two methods of representing aerosol microphysics that we employed here. The bulk aerosol treatment (referred to as One-Moment Aerosol, OMA, in Bauer et al., 2020, called “GISS-OMA” in this study) involves specification of an aerosol dry radius, and the aerosols grow hygroscopically as a function of the relative humidity (Koch et al., 2006). The more complex aerosol treatment (referred to as Multiconfiguration Aerosol TRacker of mIXing state, MATRIX, in Bauer et al., 2008, called “GISS-MATRIX” in this study) involves computation using quadrature of moments and can represent aerosol microphysical growth via condensation and coagulation. While not strictly a modal treatment, as it computes aerosol size using a quadrature of moments approach (usually termed a two-moment approach), it is most comparable to the modal treatments used in the other two models. For anthropogenic, tropospheric aerosols, a comparison between the two aerosol treatments in GISS is presented in Bauer et al. (2020). Both aerosol methods have representations of heterogeneous halogen chemistry on the aerosol surfaces and are coupled to the same chemistry scheme, which in the stratosphere includes chlorine and bromine chemistry as well as a treatment of polar stratospheric clouds (Shindell et al., 2006).

2.3 UKESM1.0

The UKESM1 Earth system model (Sellar et al., 2019) was developed jointly by the UK's Met Office and Natural Environment Research Council and consists of the physical atmosphere–land–ocean–sea ice model HadGEM3-GC3.1 (Kuhlbrodt et al., 2018) coupled to components which deal with terrestrial and ocean biogeochemistry (Wiltshire et al., 2021; Yool et al., 2021; Walters et al., 2019). Its comprehensive stratospheric–tropospheric chemistry (StratTrop) scheme is described in Archibald et al. (2020) and includes detailed tropospheric chemistry and a somewhat simplified stratospheric chemistry with “lumped” treatment of long-lived halogenated ozone-depleting substances. The model is further coupled to a modal aerosol microphysics scheme described in Mulcahy et al. (2018) and Mann et al. (2010) (GLOMAP-mode). The horizontal resolution is 1.875^∘ longitude × 1.75^∘ latitude, with 85 vertical levels up to ∼84 km on terrain-following hybrid height coordinates. An evaluation of simulations of past volcanic eruptions has been recently performed in Dhomse et al. (2020) for on older version of this model using a similar version of the aerosol scheme, as well as in Aubry et al. (2021).

2.4 Differences between the aerosol microphysics schemes used

None of the models used in this study employ sectional aerosol microphysics, whereby the size domain is divided into intervals, or bins, and the evolution of the number concentrations in each size bin is calculated separately for each species. Instead, these models employ a more simplified approach, either modal (or comparable to modal in the case of GISS-MATRIX), whereby the aerosol population is described by a number of lognormal distributions, or bulk, whereby the size distribution is prescribed; we will describe both briefly below.

CESM2 and UKESM1 use a very similar modal approach, whereby the aerosol population is described by at least three main modes (at least for sulfate) called nucleation (only for UKESM1), Aitken, accumulation and coarse, whose distribution is assumed to be lognormal with a fixed geometric standard deviation σ_g and a geometric mean diameter that can vary in a certain predefined size range: particle number and mass are transferred to the larger mode when the diameter exceeds the upper limit for that mode. In each mode, all aerosol species are considered internally mixed: that is, they are described by a single size distribution; this has been shown to lead to changes in upper-tropospheric aerosol concentrations when large quantities of stratospheric sulfate settle down, unrealistically reducing the size distribution and thus the settling velocities of species that should not interact with the sulfate (Visioni et al., 2022).

The GISS in its bulk version incorporates a much simpler aerosol parameterization, whereby the size distribution is specified for each aerosol species. There is no calculated number concentration and water uptake effects are solely dependent on relative humidity for gravitational settling purposes. GISS-MATRIX sits in between the complexity of the other two models and GISS-OMA: the sulfate population is described by two lognormal distributions (Aitken and accumulation) of fixed σ_g, but separated from other aerosol species, defined as particles with a diameter smaller or larger than 0.1 µm, respectively. Instead of tracking the size distribution itself, however, only key moments of the aerosol population such as number and mass are tracked based on the quadrature method of moments framework described in McGraw (1997), leading to a set of equations to be solved for each population as described in detail in Bauer et al. (2008). To prevent the mean diameter of the Aitken mode from approaching that of the accumulation mode too often, a transfer function that completely moves all particles from the Aitken to the accumulation mode when the two sizes approach each other in value is used (see Sect. 2.8 in Bauer et al., 2008). For a straightforward comparison, all information for the aerosol populations in each model is described in Table 1. For GISS-OMA, the only prescribed value is the bulk diameter (which for simplicity we included in the accumulation row in the table) with no assumptions regarding the geometric standard deviation; the actual value of 0.3 µm used in the current version of GISS used here is different from the one given in Koch et al. (2006) (0.4 µm).

We note that in the GISS with bulk treatment, there is one specified aerosol dry radius for all sulfate aerosols, both tropospheric and stratospheric. We used the default aerosol size, which is calibrated to represent tropospheric sulfate and the background sulfate layer in the stratosphere but is far too small for the aerosols that would result from a high stratospheric loading of sulfate. Past experiments using an earlier version of this model increased the aerosol size, resulting in a better match to stratospheric sulfate aerosols compared to observations but was far larger than should have occurred for tropospheric aerosols. This had non-negligible effects on radiative forcing, tropospheric chemistry and aerosol deposition (e.g., Kravitz et al., 2009; Pitari et al., 2014), limiting the ability to compare the SAI run with a corresponding baseline run. The approach chosen in the present study avoids these issues, but in doing so limits the applicability of the GISS-OMA simulations to SAI. Nevertheless, these simulations serve as a useful point of comparison and reveal understanding.

As shown in Table 1, the modal aerosol schemes present differences in the definition of the size ranges per each mode, and we will show that this significantly changes which modes the aerosols find themselves in. While the width of the smallest mode is similar between all three models, it is significantly different for the coarse mode, being 1.2 for CESM2 and 2 for UKESM. Considering that σ_g is a parameter that remains constant in the whole atmosphere, this difference may have arisen if the two models considered different types of aerosol size distributions to which they desired their mode to be more similar. For instance, if one considers the stratospheric size distribution of aerosols post-Pinatubo, σ_g has values lower than 1.5 (Post et al., 1997), whereas for background conditions or small eruptions, the value is closer to 1.8 (Wurl et al., 2010; Doeringer et al., 2012). As detailed in Mills et al. (2017), CESM modified its coarse mode in MAM4 to better reproduce the Pinatubo plume.

Table 1 Values for the size range (diameter, D, in µm) and geometric standard deviation σ_g for sulfate of each of the modes (when applicable) for the four participating models. Information collected from Koch et al. (2006) for GISS-OMA, Bauer et al. (2008) for GISS-MATRIX, Walters et al. (2019) for UKESM and Liu et al. (2012) with integrations from Mills et al. (2016) for CESM2. All diameters are considered here to be dry, as models calculate water uptake afterward depending on local humidity.

n/a: not applicable.

3 Experimental protocol

3.1 Baseline emission scenario

For these experiments, we selected the baseline Shared Socioeconomic Pathway (SSP) 2-4.5 (Meinshausen et al., 2020). Each experiment is run for 10 years starting on 1 January 2035. The choice of a particular background scenario is of secondary importance here, as all comparisons of the SAI responses will be performed relative to the control SSP2-4.5 simulation, thereby evaluating climate changes due solely to the increased sulfate aerosol concentrations. However, in planning for more comprehensive experiments in the future (as those detailed in Richter et al., 2022, and MacMartin et al., 2022), we will use a similar background greenhouse gas scenario.

3.2 Specifics for the SO₂ injection

In the experiments discussed in Tilmes et al. (2017), the altitudes of SO₂ injections were defined in terms of a fixed height above the annual mean tropopause. This led to two different injections altitudes, one for 15^∘ N and 15^∘ S injections at 25.0 km and one for injections at 30^∘ N and 30^∘ S at 22.4 km. Here we define a single injection altitude for all latitudes at 22 km. We chose this injection altitude for two separate reasons. First, it is closer to the upper limit of altitudes achievable by traditional aircraft (Smith, 2020) while still being sufficiently far above the tropopause for aerosols not to be removed too quickly (as opposed to 25 km, which is likely too high for practical considerations). Second, we aimed to inject SO₂ in all models at the same latitude but also in just one grid box; however, there is a challenge in prescribing a fixed injection altitude in kilometers but having a hybrid height vertical coordinate (UKESM1) or a hybrid sigma–pressure coordinate (CESM2, GISS). After performing some tests, the specific choice of 22 km ensured that all three models could always be reasonably expected to inject in the same grid box, as that altitude is not close to the vertical edges of a grid box in any of the models. In CESM2, this meant injecting at the grid box bounded by the 47 and 39 hPa pressure interfaces (with a midpoint of 43 hPa), which has an average geometric height of 21.6 km at all considered latitudes of injection and is 1.2 km thick. In UKESM1, this meant injecting at the level with a midpoint height of 21.8 km, which is roughly 1 km thick. In GISS, this meant injecting at the grid box bounded by the 43 and 31 hPa interfaces, with an approximate midpoint of 37 hPa (23.1 km assuming a 7 km scale height) and approximately 2.3 km thick.

SO₂ is injected continuously throughout the year in the same quantity in each experiment: each simulation includes the injection in only one location, either 30^∘ N, 15^∘ N, 0^∘ N, 15^∘ S or 30^∘ S, and at one longitude (180^∘ E for GISS and CESM2, 0^∘ E for UKESM1) at the altitude described above. Simulations are performed for 10 years each. The selected SO₂ injection rate in each model is 12 Tg SO₂ per year to allow for an easier comparison with past simulations with CESM1-WACCM in Tilmes et al. (2017).

4 Results

Changes in stratospheric aerosols

Figure 1 shows the evolution of the global mean changes in stratospheric aerosol optical depth (AOD, calculated for all models at 550 nm) and surface temperatures. Relatively good agreement is found between CESM2, GISS-MATRIX and UKESM1 in terms of the global mean AOD responses, whereas the global mean temperature response is similar for CESM and UKESM1, as well as for the two GISS versions separately. GISS-OMA presents AOD values larger than the other two models and GISS-MATRIX, with differences ranging from values that are 2 times larger (for equatorial and 15^∘ N injections) to 33 % larger (for 30^∘ N injections). Larger values of sulfate burden can result in a number of nonlinear effects (Niemeier and Timmreck, 2015), for instance on atmospheric dynamics and thus the latitudinal aerosol distribution (Visioni et al., 2020) or on surface climate (Simpson et al., 2019; Jiang et al., 2019). The uncertainties in the SAI process that leads from SO₂ injection to surface impacts can be separated into three different parts: (i) at each latitude how much SO₂ is needed to achieve a certain optical depth (i.e., the efficiency of SO₂ to H₂SO₄ conversion and of the removal processes), (ii) the resulting distribution of AOD under specific injection location(s) (largely driven by large-scale dynamics and mixing), and (iii) the impacts of a specific aerosol distribution on climate. Simulations with fixed SO₂ injections at fixed locations, as used in this paper, allow us to explore better points (i) and (ii). In contrast, simulations like those described in Kravitz et al. (2017), wherein the amount of SO₂ injected is adjusted each year in order to achieve some specified surface temperature goals, allow better understanding of the response to some specific pattern of AOD forcing (as the injection rate can be adjusted to achieve a desired AOD). In the following analyses we separate and discuss both the overall simulated zonal mean response and that normalized by the global magnitude of the response to highlight these different contributors to the overall uncertainty.

Figure 1 Time series of global mean monthly mean changes in stratospheric AOD (a–e) and surface temperature (f–j) resulting from the single-point injections at 30^∘ N (a, f), 15^∘ N (b, g), 0^∘ N (c, h), 15^∘ S (d, i) and 30^∘ S (e, j) compared to the corresponding baseline SSP2-4.5 simulation in each of the models. The dashed line in all panels represents the month after which averages are performed for the following figures (i.e., the last 7 years of simulations).

Figure 2 shows the latitudinal distributions of the simulated stratospheric AOD. The inter-model spread is different depending on the injection location in terms of both overall magnitude and the spatial distribution. Equatorial injections show by far the greatest differences, with UKESM1 simulating twice as much AOD in the tropics compared to CESM2 and the two GISS versions differing by up to a factor of 5 between them. Similar differences between UKESM1 and CESM2 have also been documented for the GeoMIP G6 experiment in Jones et al. (2021). Some of these differences are driven by differences in the global mean values of AOD: as mentioned above, AOD in GISS-OMA is over 2 times larger than in GISS-MATRIX for the equatorial injection case. In the right panels of Fig. 2 the zonal mean AOD values are scaled by the respective global mean values; this highlights the differences in the latitudinal distribution of the responses. This way, the difference in the normalized magnitudes of the equatorial peak (defined here as the average between 5^∘ N and 5^∘ S) between CESM2 and UKESM1 is reduced from 1.8 to 1.4; similarly, the two GISS versions show much more similar results, especially in the tropical region. In general, injections in the Southern Hemisphere (SH) show larger inter-model differences compared to the Northern Hemisphere (NH) injections, especially at high latitudes. In the 30^∘ S case, global mean values are more than twice as large in GISS-OMA (0.25) compared to UKESM1 (0.12), and the SH high-latitudinal AOD differs by a factor of 3. In the NH, the largest difference is 60 % between UKESM1 (0.12) and GISS-OMA (0.20) in the 30^∘ N case, and high-latitudinal AOD differs by a factor of 2. When normalized, the differences are largely reduced, but a clear differentiation remains for UKESM1, which shows a larger tropical confinement of aerosols than the other models. There are also some notable differences in the interannual variability between models: while both GISS versions and CESM2 have standard deviations that range between 1 % and 5 % of the mean values (depending on the precise latitude), UKESM1 shows a much larger variability (up to 33 % in the 30^∘ N injection case at 60–90^∘ N). In future longer simulations this observation would need to be considered when also considering the interannual variability of the surface changes.

The differences in AOD can be better understood by looking at the changes in aerosol mass mixing ratios (Fig. 3) and the overall sulfate column burden (Fig. 4). Compared to the other models, UKESM1 shows a much stronger confinement of the aerosol mass in the tropical pipe (for equatorial injections, the peak is 1.5 µg SO₄ kg⁻¹ of air, whereas it is 1.1 µg SO₄ kg⁻¹ of air for CESM2 and GISS-OMA and 0.5 for GISS-MATRIX) and far less transport in the opposite hemisphere for both 15^∘ injections. On the other hand, CESM2 shows a much stronger poleward transport of sulfate aerosols, as indicated by aerosol values twice as large for 30^∘ injections compared to other models considering values above 60^∘ in latitude. All of these characteristics can largely be explained by the characteristics of the climatological transport in each of these models, as discussed in depth in Part 2.

GISS-OMA also shows a larger poleward transport in all cases, but the underlying cause for this might be different between GISS and CESM2. As analyzed in depth in Part 2, the baseline stratospheric dynamics are very similar between the two GISS realizations, and hence the differences in the overall aerosol distribution cannot be found there. Analyses of the stratospheric aerosol surface area density (SAD, shown in Part 2 due to its importance for stratospheric ozone chemistry) indicates a far larger number of particles in GISS compared to the other models. The SAD in GISS-OMA is 3 times as large poleward of 60^∘ compared to CESM2, even for mass values that only differ by 30 % (for instance, in the 30^∘ N injection case). GISS-OMA simulates much smaller particles (the prescribed bulk dry radius is 0.15 µm) that imply higher residence times, as gravitational settling is much lower (Visioni et al., 2022), so the particles are more easily transported towards the pole even in the presence of a similar dynamical regime compared with the other GISS realization. The presence of aerosols at higher altitudes even close to the poles (where the large-scale circulation should be transporting aerosols downward) further supports this observation. Smaller aerosols are also more efficient scatterers, explaining why the optical depth differences are larger than the mass differences.

Finally, panels (f)–(h) in Fig. 4 give an overview for all cases of the ratio between the overall mass of the produced aerosols (shown alone in Fig. 4f) and the injected amount of SO₂. This ratio, shown in Fig. 4g, represents the lifetime of the added sulfate: as we are in a steady state, wherein no new mass is added to the global burden, this lifetime can be calculated as the burden divided by its constant source (the injection) (Visioni et al., 2018). All models show larger amounts of mass for injections closer to the Equator due to the tropical confinement increasing the aerosol lifetime. This effect is less noticeable for GISS-OMA, which shows far less confinement. However, models show a distinct difference in the ratio between the overall mass of SO₄ and resulting AOD (Fig. 4h): this indicates substantial inter-model differences in the average size of the aerosol, with obvious consequences for the efficiency of the cooling that we will analyze next.

A comparison of the effective radius (R_eff) is shown in Fig. 5. The values for R_eff are indirectly derived for all models as R_eff is not a direct output for either UKESM1 or GISS. Therefore, we use the common available output of mass mixing ratio (χ) and number concentration (N) for each mode in each model to derive the mean radius r_i. We calculate the mean volume v_i using

$$\begin{array}{}\text{(1)}& v_i={\displaystyle \frac{{\mathit{\chi }}_i}{{\mathit{\rho }}_{\mathrm{sulfate}}\cdot N_i}},\end{array}$$

where ρ_sulfate is the density of sulfate as considered in each model in kilograms per cubic meter (kg m⁻³) (usually 1770 kg m⁻³ as in Liu et al., 2012). The mean radius is then derived considering that the aerosol population is assumed to be a lognormal distribution with a geometric standard deviation σ_g and is thus connected to the mean volume by

$$\begin{array}{}\text{(2)}& v_i={\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\cdot \mathit{\pi }\cdot r_i^{\mathrm{3}}\cdot \exp \left({\displaystyle \frac{\mathrm{9}}{\mathrm{2}}}\cdot {\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_{\mathrm{g}}\right)\right).\end{array}$$

Finally, R_eff is calculated considering the definition

$$\begin{array}{}\text{(3)}& R_{\mathrm{eff}}={\displaystyle \frac{{\sum }_i{r}_i^{\mathrm{3}}{N}_i}{{\sum }_i{r}_i^{\mathrm{2}}{N}_i}},\end{array}$$

where i is the number of modes considered and is therefore 3 for CESM2 (all modes relevant to sulfate in MAM4, Liu et al., 2016), 4 for UKESM1 and 2 for GISS-MATRIX. The single wet radii for each mode are shown in Fig. S1 in the Supplement. The validity of our derivation is ensured by the comparison of the derived values of R_eff in CESM2 and the ones obtained as a direct output from the CESM2 simulations (not shown); while the off-line derivation leads to slightly smaller R_eff values than calculated online in the model, the overall results are similar enough for a confident comparison. A number of important features become apparent: the population-weighted radius in GISS-MATRIX is almost always larger than the GISS-OMA one, confirming our observation of the differences in high-latitudinal concentration. CESM2 also presents a radius always larger than the GISS-MATRIX and UKESM1 one, which explains most of the discrepancies between AOD and column mass between the CESM2 and GISS models: even if CESM2 shows larger mass concentrations, smaller particles such as those in GISS are more efficient scatterers; thus, GISS AOD is more comparable to CESM2. If dynamics were similar between the two models, one would also expect less mass in CESM2 compared to GISS due to reduced lifetime caused by the increased gravitational settling, but in this case, the dynamical differences dominate in determining the final mass concentrations; these differences are discussed in more depth in Part 2 and indicate a substantially stronger climatological shallow branch of the Brewer–Dobson circulation in CESM2 compared to other models (see Fig. 2 in Part 2). UKESM1 has intermediate radii between CESM2 and GISS-MATRIX.

Looking at the mean radii for each mode (Fig. S1), differences in the microphysical approaches become more evident between the models, also based on the values given in Table 1. Given that UKESM1 defines the coarse mode as particles larger than 0.5 µm, all of the stratospheric aerosols are found in the accumulation mode (as they are in GISS), whereas CESM2, which defines the coarse-mode threshold at around 0.4 µm, simulates most of the sulfate aerosols in that mode and very few in the accumulation mode. This particular observation should not significantly influence the aerosols' behavior in the stratosphere, since the coagulation and condensation processes happen in all modes. It could, however, have some impacts on the tropospheric interactions of sulfate with cloud nuclei; for instance, some models (such as CESM2) treat coarse-mode and accumulation-mode aerosols differently when calculating the presence of heterogeneous and homogeneous ice nuclei (see, for example, Visioni et al., 2022, for a discussion of undesired side effects this might have on ice nucleation rates in the upper troposphere in CESM1(WACCM)).

Figure 2 (a–e) Zonal and annual mean increase in stratospheric AOD resulting from single-point injections at 30^∘ N (a), 15^∘ N (b), 0^∘ N (c), 15^∘ S (d) and 30^∘ S (e), averaged over the last 7 years of simulation. (f–j) Zonal mean increase normalized to the global value of AOD in each simulation in the respective experiment; shading of the same color represents the internal variability (1σ standard deviation) over the same period of time. For each experiment, the global mean value used in the normalization is reported between the panels on the left and right. Dashed vertical lines indicate the locations of the SO₂ injections.

Figure 3 Zonal and annual mean increase in SO₄ mass mixing ratio (in 10⁻⁷ kg SO₄ kg⁻¹ of air) for CESM2 (first column), UKESM1 (second column), GISS-MATRIX (third column) and GISS-OMA (fourth column) as well as all injection locations at 30^∘ S (first row), 15^∘ S (second row), 0^∘ N (third row), 15^∘ N (fourth row) and 30^∘ N (fifth row). Hatched areas indicate locations where changes are not statistically significant (using a two-sided t test at the 95 % confidence level).

Figure 4 (a–e) Zonal and annual mean increase in total SO₄ column burden for each experiment, 30^∘ N (first row), 15^∘ N (second row), 0^∘ N (third row), 15^∘ S (fourth row) and 30^∘ S (fifth row). (f) Total global increase in stratospheric sulfate burden (in Tg S). (g) Stratospheric sulfate lifetime (in months). (h) Global increase in column sulfate burden, normalized by the obtained global mean increase in stratospheric AOD.

Figure 5 Effective radius (in 10⁻⁶ m) in CESM2 (left), UKESM1 (center) and GISS-MATRIX (right) simulations for each experiment, 30^∘ S (first row), 15^∘ S (second row), 0^∘ N (third row), 15^∘ N (fourth row) and 30^∘ N (fifth row).

The resulting zonal mean cooling achieved by the AOD increase in each simulation is shown in Fig. 6, including both the absolute temperature changes and values normalized by the global means. Differences between models are still evident, but a direct correlation with the AOD changes in Fig. 2 is not present; for instance, the general temperature responses between the two GISS realizations are much more similar than their simulated AOD, as was already visible in Fig. 1. This is because temperature responses compound many sources of uncertainty, including those related not only to the magnitude and spatial pattern of the simulated stratospheric AOD (thus related to the model microphysics and transport), but also to the resulting radiative, dynamical and chemical processes operating in the troposphere and lower stratosphere (including the associated changes in ozone and water vapor, which are discussed in Part 2) that modulate how much cooling AOD produces and how the regional patterns of cooling differ. We find very large differences between the models in the magnitude of the global mean cooling per unit AOD (Fig. 7). CESM2 and UKESM1 show, on average, 4.7 and 5.4 K cooling per unit of AOD, respectively. This is roughly in line with the multi-model average of the GeoMIP G6sulfur experiment (six models, resulting in 4.6 K per unit AOD for injections between 10^∘ N and 10^∘ S in GeoMIP – Visioni et al., 2021 – compared to the 3.9 and 4.8 K cooling per unit AOD in CESM2 and UKESM1, respectively, for the equatorial injections here).

In contrast, the GISS models show significantly larger cooling per unit AOD, with 8.6 K on average in GISS-OMA and 16.7 K on average in GISS-MATRIX; the latter is more than 3 times as high as in CESM2. As noted for Fig. 1, this large normalized difference is mainly due to differences in the global cooling produced between the two GISS realizations (the sensitivity to the aerosol amount); even if there is a much smaller amount of AOD simulated in GISS-MATRIX, it results in a similar level of surface cooling as is simulated in GISS-OMA under much larger AOD values. The cause for the discrepancies between the two GISS realizations is most likely found in similar tuning procedures for the model (described for GISS-E2.1 in Kelley et al., 2020); while these mainly target the global radiative balance, they often also try to constrain the aerosol forcing (this is described, for instance, in Schmidt et al., 2017, where both GISS and CESM1 tuning procedures are explained in some detail). The two GISS versions used here show large differences in the baseline aerosol optical depth from sulfate: in the midlatitudinal NH, where anthropogenic sulfate aerosol production is maximized, GISS-MATRIX shows an average AOD of 0.03, just like CESM2, while GISS-OMA shows an average AOD of 0.11 (see Fig. S2). Further investigations will be warranted in the baseline GISS simulations to understand if this difference is relevant for the observed differences in the sensitivity to stratospheric aerosols between the two versions.

Another possible contributor to the differences in the magnitude of global cooling simulated in GISS could be the substantial tropospheric ozone reduction simulated in the two model versions (see Sect. 3.3.4 of Part 2). As a greenhouse gas, reduction in tropospheric ozone would thus act to magnify the cooling from the direct aerosol forcing, and this effect is generally not present in CESM2 or UKESM1.

The normalized values in the right panels of Fig. 6 show some level of consistency in the temperature response amongst the models. In particular, the models agree that equatorial injections tend to cool the NH more than the SH and that the cooling produced by the NH injections is strongest in the NH (this is particularly marked in CESM2 for the 15^∘ N injection), whereas the cooling produced by the SH injections is more evenly distributed between the two hemispheres. In contrast, the models tend to disagree with respect to responses in the NH high latitudes much more than for the SH high latitudes; this can be explained by the role of dynamical atmospheric variability and sea ice variability, which in the NH both may contribute to increasing model spread (Banerjee et al., 2021).

Figure 6 (a–e) Zonal and annual mean change in surface temperatures (K) resulting from the single-point injections at either 30^∘ N (a), 15^∘ N (b), 0^∘ N (c), 15^∘ S (d) or 30^∘ S (e), averaged over the last 7 years of simulation. (f–j) Zonal mean changes normalized with the corresponding global mean values in each experiment; shading of the same color represents the interannual variability (1σ) over the same period of time. For each experiment, the global mean value used in the normalization is reported between the panels on the left and right.

Figure 7 Global mean cooling achieved in each experiment and injection location per unit increase in the associated global mean stratospheric AOD.

Finally, in Fig. 8 we show the zonal mean changes in precipitation. In this case, we show both the magnitude of the change in millimeters per day (mm d⁻¹) and the corresponding percent changes with respect to the baseline precipitation pattern in each model. The precipitation response also depends on the climatological precipitation, which differs more between the models than the climatological temperature does (not shown). We find general qualitative agreement between the models in the tropics. All models show shifts in tropical precipitation depending on the hemisphere of injection, with decreased precipitation in that hemisphere and increased precipitation in the opposite one. The equatorial injections in turn result in an overall reduction of tropical precipitation. Changes at middle and high latitudes, which are much smaller in absolute terms than those in the tropics, can be as high in terms of percentage changes as those in the tropics, but the models disagree on the sign and magnitude of the responses there far more than they do in the tropics. We also show percent in order to highlight the significance of the absolute changes based on the overall received precipitation. This, for instance, helps to show that, compared to the baseline values, the very small changes observed at high latitudes may also be significant.

In general, precipitation changes depend on even more factors than temperature changes (Kravitz et al., 2013; Tilmes et al., 2013). Differences in the inter-hemispheric cooling patterns produce shifts in the precipitation by moving the Intertropical Convergence Zone (ITCZ, see below and, e.g., Ridley et al., 2015); this is compounded by differences in the cloud responses (Smyth et al., 2017) and by the large-scale dynamical changes resulting from the aerosol-induced stratospheric heating (see Sect. 3.2 to 3.5 in Part 2; also Simpson et al., 2019). Furthermore, precipitation responses tend to be regionally dependent: for instance, previous G6 GeoMIP simulations showed consensus in the winter precipitation response over Europe (with increased precipitation in northern Europe and decreased over southern Europe associated with the SAI-induced positive phase of the North Atlantic Oscillation); however, the same models showed little consensus over the North American continent (Jones et al., 2022).

Figure 8 (a–e) Zonal and annual mean changes in precipitation (mm d⁻¹) resulting from the single-point injections at either 30^∘ N (a), 15^∘ N (b), 0^∘ N (c), 15^∘ S (d) or 30^∘ S (e), averaged over the last 7 years of simulation. Panels (f)–(j) are as (a)–(e), but with percent changes obtained by normalizing to the baseline zonal mean precipitation in each model.

Figure 9 shows changes in the position of ITCZ, which we define here as the latitude near the Equator where the meridional mass streamfunction at 500 hPa changes sign. The models qualitatively agree with respect to the simulated ITCZ shifts, despite the significant differences in the climatological ITCZ locations (shown on the left side of Fig. 9). The larger simulated ITCZ shifts for the NH injections compared to the SH injections are also consistent with the different cooling patterns, with stronger and more NH-confined cooling for the 30 and 15^∘ N injections. A similar behavior of the ITCZ was also previously found for volcanic eruptions occurring in the equatorial or extra-equatorial regions (Trenberth and Dai, 2007; Oman et al., 2006; Haywood et al., 2013; Ridley et al., 2015) as well as for previous single-point SAI simulations (Haywood et al., 2013).

Figure 9 Annual mean changes in the latitude of the ITCZ, defined here as the latitude near the Equator where the meridional mass streamfunction at 500 hPa changes sign (defined as the midpoint between the two grid boxes). Results are shown for each injection location (30^∘ S, 15^∘ S, 0^∘ N, 5^∘ N and 30^∘ N) and each model. Error bars indicate the σ standard deviation for the last 7 years of the SAI simulation. Also shown on the left-hand side is the climatological ITCZ position in the corresponding SSP2-4.5 simulation.

5 From single-point injections to multi-target simulations: designing the feedback algorithm

One of the goals of this exercise is to reproduce in multiple climate models a geoengineering strategy similar to that described in Kravitz et al. (2017), i.e., wherein a feedback algorithm determines the amount of SO₂ to be injected each year at each of the specified locations in order to achieve a set of predetermined surface temperature targets. As discussed in Kravitz et al. (2016) and Zhang et al. (2022), year-round injections at the four injection locations described here (30^∘ N, 15^∘ N, 15^∘ S and 30^∘ S) may allow for the management of three independent degrees of freedom in surface temperature (global mean temperature, temperature difference between the two hemisphere, temperature difference between the pole and the tropics). Combining varying rates of injections at these four locations produces distinct patterns of AOD.

One method of describing and decomposing these AOD patterns into simpler ones is to project the zonal mean AOD onto the first three Legendre polynomials evaluated as a function of sine of latitude (hereafter, ψ), which equate to three idealized patterns of optical depth achievable. Here ℓ₀=1 indicates a pattern of AOD uniform at all latitudes, which therefore produces a uniform cooling; ℓ₁=sin (ψ) indicates a pattern of AOD that is greater in one hemisphere and smaller in another, which can therefore cool one hemisphere more than another; and ${\mathrm{\ell }}_{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{3}{\sin }^{\mathrm{2}}\right(\mathit{\psi })-\mathrm{1})$ indicates a quadratic pattern that is smaller close to the Equator and greater at high latitude, therefore cooling those latitudes more. MacMartin et al. (2017) demonstrated that these three patterns are achievable through the combination of injections at different locations: injecting at the same time at 15^∘ N and 15^∘ S produces ℓ₀; injecting at 15 and 30^∘ N (or S) produces ℓ₀+ℓ₁ (${\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}$, more AOD in the Northern Hemisphere) or ℓ₀−ℓ₁ (${\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}$, more AOD in the Southern Hemisphere), respectively; and injecting at 30^∘ N and 30^∘ S produces ℓ₀+ℓ₂. Therefore, injections at these four latitudes can be combined (assuming linearity) to achieve desired combinations of ℓ₀, ℓ₁ and ℓ₂, subject to the constraints ℓ₀≥0, ℓ₂≥0 and ${\mathrm{\ell }}_{\mathrm{0}}\ge \left|{\mathrm{\ell }}_{\mathrm{1}}\right|+{\mathrm{\ell }}_{\mathrm{2}}$. These constraints derive from the fact that a negative injection (and a negative pattern of AOD) is not possible.

This truncated Legendre decomposition of the zonal mean AOD into ℓ₀, ℓ₁ and ℓ₂ is not the only way to numerically define an SAI strategy, but it has been found to be sufficient to independently control the three degrees of freedom in surface temperature: the global mean temperature, $T_{\mathrm{0}}=\frac{\mathrm{1}}{A}{\int }_{\mathit{\psi }}T\left(\mathit{\psi }\right)\cos \left(\mathit{\psi }\right)\mathrm{d}\mathit{\psi }$ (where T(ψ) is the zonal mean surface temperature and $A={\int }_{\mathit{\psi }}\cos \left(\mathit{\psi }\right)\mathrm{d}\mathit{\psi }$ is a normalization factor), controlled by ℓ₀; the pole-to-pole temperature gradient $T_{\mathrm{1}}=\frac{\mathrm{1}}{A}{\int }_{\mathit{\psi }}\sin \left(\mathit{\psi }\right)T\left(\mathit{\psi }\right)\cos \left(\mathit{\psi }\right)\mathrm{d}\mathit{\psi }$, controlled by ℓ₁; and the Equator-to-pole temperature gradient $T_{\mathrm{2}}=\frac{\mathrm{1}}{A}{\int }_{\mathit{\psi }}\frac{\mathrm{3}{\sin }^{\mathrm{2}}\left(\mathit{\psi }\right)-\mathrm{1}}{\mathrm{2}}T\left(\mathit{\psi }\right)\cos \left(\mathit{\psi }\right)\mathrm{d}\mathit{\psi }$, controlled by ℓ₂. Similarly, the definition of the four injection latitudes described here is not exhaustive in the sense that it is not the only way to produce the three described aerosol patterns, as there could be other latitudes of injection that once combined could produce ℓ₀, ℓ₁ and ℓ₂ with smaller residuals. Our choice is thus only intended to illustrate the optimization method.

In this section, we present the calculations necessary for this strategy to be implemented in the three climate models considered in this study; similar analyses for CESM1(WACCM) are performed and explained in detail in MacMartin et al. (2017). In order to design such a feedback control algorithm, it is first necessary to quantify the relationships between injections at each of the four latitudes (q=[q_30 S q_15 S q_15 N q_30 N]^T, where q is the amount of injection per year in Tg SO₂) as well as the patterns of AOD they create (ℓ=[ℓ₀ ${\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}$ ${\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}$ ℓ₂]^T) and then the relationship between AOD (ℓ) and temperature (T₀, T₁, and T₂). The first relationship can be defined by two matrices M and F such that $q={\mathbf{MF}}^{-\mathrm{1}}\mathrm{\ell }$. For each of the four models, using the values shown in Fig. 2, this relationship is computed to be as follows.

$$\begin{array}{}\text{(4)}& \begin{array}{r}\text{CESM: }\left[\begin{array}{c}{q}_{\mathrm{30}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{N}}\\ q_{\mathrm{30}\mathrm{N}}\end{array}\right]=\left[\begin{array}{c}\mathrm{25}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+\mathrm{38}{\mathrm{\ell }}_{\mathrm{2}}\\ \mathrm{39}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{45}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}\\ \mathrm{34}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{59}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}\\ \mathrm{15}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+\mathrm{43}{\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\\ \text{UKESM: }\left[\begin{array}{c}{q}_{\mathrm{30}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{N}}\\ q_{\mathrm{30}\mathrm{N}}\end{array}\right]=\left[\begin{array}{c}\mathrm{56}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+\mathrm{50}{\mathrm{\ell }}_{\mathrm{2}}\\ \mathrm{37}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{22}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}\\ \mathrm{36}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{18}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}\\ \mathrm{49}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+\mathrm{46}{\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\\ {\mathrm{GISS}}_{\mathrm{OMA}}:\left[\begin{array}{c}{q}_{\mathrm{30}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{N}}\\ q_{\mathrm{30}\mathrm{N}}\end{array}\right]=\left[\begin{array}{c}\mathrm{6}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+\mathrm{26}{\mathrm{\ell }}_{\mathrm{2}}\\ \mathrm{23}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{35}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}\\ \mathrm{20}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{44}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}\\ \mathrm{0}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+\mathrm{29}{\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\\ {\mathrm{GISS}}_{\mathrm{MATRIX}}:\left[\begin{array}{c}{q}_{\mathrm{30}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{S}}\\ q_{\mathrm{15}\mathrm{N}}\\ q_{\mathrm{30}\mathrm{N}}\end{array}\right]=\left[\begin{array}{c}\mathrm{38}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+\mathrm{46}{\mathrm{\ell }}_{\mathrm{2}}\\ \mathrm{44}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{47}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}\\ \mathrm{39}({\mathrm{\ell }}_{\mathrm{0}}-({\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}+{\mathrm{\ell }}_{\mathrm{2}}\left)\right)+\mathrm{69}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}\\ \mathrm{22}{\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}+\mathrm{49}{\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\end{array}\end{array}$$

These values can be obtained by separately considering the injection amounts necessary to obtain the four AOD patterns described above and solving the linear least-squares problem considering as a constraint that no negative injection or AOD values are possible. This is shown in Fig. 10 for all models, with the dashed black lines being the idealized patterns; the linear combination of the coupled injection locations in each of the four cases, which results in the colored lines in each plot, successfully manages to produce a pattern similar to the idealized one, indicating that by injecting in the locations presented in the tables below each panel the quantities indicated, those patterns are achievable. Models all agree on the need to deploy quasi-symmetrical amounts to obtain the ℓ₀ and ℓ₀+ℓ₂ patterns (with at most a 20 % difference in some cases), whereas they differ largely on how to obtain the ℓ₀+ℓ₁ and ℓ₀−ℓ₁ patterns. For instance, to obtain ℓ₀+ℓ₁, CESM2 requires 4 times as much injection at 15^∘ N than it does at 30^∘ N, whereas UKESM1 requires 2 times as much injection at 30^∘ N than it does at 15^∘ N, and GISS requires no injection at 30^∘ N at all. The residuals (the difference between the solution of the linearized system and the idealized pattern) also show that how close the models get to the idealized patterns differs, with CESM2 being twice as low as any other models in all cases. It is important to note that the results obtained with the current version of CESM2 (CESM2-WACCM6) are very similar but not identical (<10 %) to those obtained and discussed in MacMartin et al. (2017) using CESM1(WACCM). In that case, injections were at different altitudes (25 km for 15^∘ N and S, 22 km for 30^∘ N and S), so some differences are to be expected.

Figure 10 Linear least-squares solution to obtain an ℓ₀-shaped AOD (black dashed line) using only 15^∘ N and 15^∘ S injections (a), an ℓ₀+ℓ₁ (${\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{N}}$, black dashed line, increased AOD in the Northern Hemisphere) one using only 15 and 30^∘ N injections (b), an ℓ₀−ℓ₁ (${\mathrm{\ell }}_{\mathrm{1}}^{\mathrm{S}}$, black dashed line, increased AOD in the Southern Hemisphere) one using only 15 and 30^∘ S injections (c), and an ℓ₂ (black dashed line, increased AOD at high latitudes compared to low latitudes) one using only 30 and 30^∘ S injections (d). Values at the bottom of each plot (in Tg SO₂ yr⁻¹) indicate the required injection, assuming linearity, based on the AOD distribution in Fig. 2. The residual of each solution is also included.

Using the patterns of temperature response shown in Fig. 6 combined with the patterns in AOD, we can derive the second set of relationships between the AOD projections ℓ and T₀ as well as T₁ and T₂. For the four models, this leads to the following.

$$\begin{array}{}\text{(5)}& \begin{array}{r}\text{CESM: }\left[\begin{array}{c}{T}_{\mathrm{0}}\\ \mathrm{3}{T}_{\mathrm{1}}\\ \mathrm{5}{T}_{\mathrm{2}}\end{array}\right]\simeq \left[\begin{array}{cccc}\mathbf{-}\mathbf{4.1}& \mathit{0}& \mathit{0}& \\ -\mathrm{3.0}& \mathbf{-}\mathbf{3.9}& \mathit{0}\\ -\mathrm{1.5}& -\mathrm{1.6}& \mathbf{-}\mathbf{0.5}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ell }}_{\mathrm{0}}\\ {\mathrm{\ell }}_{\mathrm{1}}\\ {\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\\ \text{UKESM: }\left[\begin{array}{c}{T}_{\mathrm{0}}\\ \mathrm{3}{T}_{\mathrm{1}}\\ \mathrm{5}{T}_{\mathrm{2}}\end{array}\right]\simeq \left[\begin{array}{cccc}\mathbf{-}\mathbf{4.6}& \mathit{0}& \mathit{0}& \\ -\mathrm{2.6}& \mathbf{-}\mathbf{3.6}& \mathit{0}\\ -\mathrm{0.2}& -\mathrm{1.5}& \mathbf{-}\mathbf{1.3}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ell }}_{\mathrm{0}}\\ {\mathrm{\ell }}_{\mathrm{1}}\\ {\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\\ {\mathrm{GISS}}_{\mathrm{OMA}}:\left[\begin{array}{c}{T}_{\mathrm{0}}\\ \mathrm{3}{T}_{\mathrm{1}}\\ \mathrm{5}{T}_{\mathrm{2}}\end{array}\right]\simeq \left[\begin{array}{cccc}\mathbf{-}\mathbf{5.6}& \mathit{0}& \mathit{0}& \\ -\mathrm{0.9}& \mathbf{-}\mathbf{4.4}& \mathit{0}\\ \mathrm{0.7}& -\mathrm{1.0}& \mathbf{-}\mathbf{4}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ell }}_{\mathrm{0}}\\ {\mathrm{\ell }}_{\mathrm{1}}\\ {\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\\ {\mathrm{GISS}}_{\mathrm{MATRIX}}:\left[\begin{array}{c}{T}_{\mathrm{0}}\\ \mathrm{3}{T}_{\mathrm{1}}\\ \mathrm{5}{T}_{\mathrm{2}}\end{array}\right]\simeq \left[\begin{array}{cccc}\mathbf{-}\mathbf{9.2}& \mathit{0}& \mathit{0}& \\ -\mathrm{4.2}& \mathbf{-}\mathbf{6.8}& \mathit{0}\\ -\mathrm{1.4}& -\mathrm{2.4}& \mathbf{-}\mathbf{6}\end{array}\right]\left[\begin{array}{c}{\mathrm{\ell }}_{\mathrm{0}}\\ {\mathrm{\ell }}_{\mathrm{1}}\\ {\mathrm{\ell }}_{\mathrm{2}}\end{array}\right]\end{array}\end{array}$$

This result indicates that the main perturbations produced on T₀, T₁ and T₂ result from the respective AOD patterns, as the larger numbers are on the main diagonal. A controller algorithm based on the values shown in Eqs. (4) and (5) can, in a strategy that aims to achieve a certain climatic target, read the annually averaged surface temperatures at the end of 1 year, determine the distance between those and the target ones, and determine how the cooling can be achieved through the three AOD patterns following Eq. (5). Then, with Eq. (4), the AOD pattern necessary can be related to four different amounts of injections at the four locations. For CESM2, these results have already been applied to produce two datasets (described in Richter et al., 2022, and MacMartin et al., 2022), and similar efforts are ongoing for UKESM and GISS.

Overall, the results shown in this section confirm that the methods applied in MacMartin et al. (2017) to produce the results in Kravitz et al. (2017) may be applicable to other climate models and that the inter-model differences in transport and microphysics can, in this kind of optimization strategy, be overcome by modifying the ratio between injection amounts at different locations. However, it is also clear that, as indicated by the differences in the residuals shown in Fig. 10, the use of the four specific latitudes proposed in MacMartin et al. (2017) might not result in the most optimal set of latitudes in all models, depending on the specifics of the stratospheric circulation. The results of UKESM, for instance, suggest that, in that model, the coordinates 15^∘ N and 15^∘ S are too close to the tropical pipe and therefore still confine the aerosols too much at low latitudes, as opposed to CESM2 and GISS. Therefore, based on UKESM results, a choice of latitudes further from the Equator (such as 20^∘) may result in a better optimization. Similarly, all models seem to suggest that the successful obtainment of ℓ₂ might be hard to achieve with only 30^∘ injections and that some higher-latitude injections might be necessary. It is thus worth reiterating that the specific choices of latitudes here should not be read as the only or optimal choice but merely as a way to start exploring the available space (see also Zhang et al., 2022, for a broader exploration of this issue).

6 Conclusions

In this work, we have shown the results of a first systematic intercomparison of climate responses to fixed single-point SO₂ injections in different climate models with interactive aerosol microphysics and comprehensive stratospheric chemistry. In particular, we used CESM2-WACCM6, UKESM1.0 and GISS-E2.1-G with bulk aerosol microphysics and GISS-E2.1-G with two-moment aerosol microphysics, as well as a set of simulations injecting a fixed quantity of SO₂ at five different latitudes (30^∘ N, 15^∘ N, 0^∘ N, 15^∘ S and 30^∘ S) and at 22 km of altitude. The same protocol was used in all models to determine similarities and differences in the resulting stratospheric sulfate aerosol distributions, alongside the resulting atmospheric and surface climate responses. Similar simulations had previously only been performed with one climate model, CESM1(WACCM) (MacMartin et al., 2017). Our multi-model simulations therefore serve a dual purpose: (i) to evaluate the responses to off-equatorial SO₂ injections in multiple climate models by understanding similarities and differences between different independent climate models and thus to isolate sources of uncertainty in model SAI responses and identify future areas of improvements and (ii) to lay the basis for a future intercomparison between models using a feedback algorithm capable of achieving multiple temperature targets (MacMartin et al., 2017) similar to the one used in Tilmes et al. (2018a) to produce the Geoengineering Large Ensemble.

The simulated changes in stratospheric aerosol optical depth were analyzed in terms of both their absolute magnitudes and the values normalized per unit of global mean AOD achieved. The latter assumes that the latitudinal distribution of simulated AOD remains similar under varying injection rates. Previous analyses of the GLENS results showed that this assumption holds, unless very large injections are considered that might noticeably modify the stratospheric circulation and prevent the aerosols from reaching the higher latitudes (Visioni et al., 2020). Our results in this work showed that in order to achieve similar global AOD values, the models considered would require different amounts of Tg SO₂ yr⁻¹ injected, with the largest differences for the model with a different, simpler aerosol treatment (GISS-E2.1-G with bulk microphysics). By analyzing the results in terms of the normalized distributions, which bypasses the bias due to different simulated global mean AOD values under the same injection rate, the largest inter-model differences were found in the tropics (especially for the equatorial injections) and at very high latitudes, with models disagreeing on the amount of aerosol transported poleward.

Using a similar separation between the absolute magnitudes and the normalized latitudinal changes we analyzed the resulting surface temperature responses. As before, a large discrepancy between the two GISS versions and the other two models was identified for the global mean response, with the former showing a global cooling per unit AOD roughly a factor of 2 (for the bulk aerosol version) or 3 (for the modal aerosol version) larger than that simulated by either CESM2-WACCM6 or UKESM1.0. When the surface temperature responses are normalized, however, models generally show good agreement amongst them in terms of the overall latitudinal distribution of the temperature changes, and similarly they show good agreement in the response of the latitudinal changes in precipitation.

Uncertainties in the projected SAI responses closely track those reported for simulations of past explosive volcanic eruptions (Clyne et al., 2021). Our results show that a large fraction of these uncertainties arises from discrepancies in the SO₂-to-sulfate aerosol conversion and their subsequent growth in the three models with three different model aerosol schemes, although differences in the stratospheric circulation also play an important role, as shown in Part 2. In addition, comparison of the two sets of GISS simulations using either two-moment or bulk treatment of aerosol microphysics showed important differences in terms of the simulated AOD and the overall cooling produced per unit AOD (Fig. 7), as well as in terms of the resulting stratospheric response (Part 2), highlighting the importance of detailed treatment of microphysical processes. In agreement, a recent study from Laakso et al. (2022) using ECHAM-HAMMOZ showed that the choice of aerosol scheme (modal versus sectional) can lead to large discrepancies in the resulting radiative forcing due to the injection of different quantities of SO₂ (up to twice as large for some scenarios). Our GISS results here further confirm that the choice of aerosol microphysical scheme, all else being equal, can influence the overall amount of cooling, although the results are still even more influenced by the choice of climate model.

This work has shown that it would be feasible to replicate a global-scale injection strategy, such as the one used in GLENS, in multiple climate models; if the amount of SO₂ injected can be controlled year by year, models appear capable of reproducing a similar scaled surface temperature response. This thus ensures that a strategy including injections at 30^∘ N, 15^∘ N, 15^∘ S and 30^∘ S would potentially be able to maintain the three defined temperature targets (i.e., global mean surface temperature; Equator-to-pole and inter-hemispheric surface temperature gradients), albeit with differences in the injection magnitudes between models. A future comparison of the results of such an experiment using different climate models maintaining similar temperature targets would help identify different sources of uncertainty in the modeled response to SAI compared to fixed-point injection simulations, for instance regarding the behavior of the Atlantic Meridional Overturning Circulation (Tilmes et al., 2020), the North Atlantic Oscillation (Jones et al., 2022) or impacts on middle atmospheric composition, including stratospheric ozone (Tilmes et al., 2022). When it comes to climate impacts, focusing on particular large-scale temperature targets rather than analyzing the direct outcomes of fixed single-point injections shifts, in a way, the focus from one kind of uncertainty to another. If the location and amount of injection can be modified, some of the stratospheric uncertainties discussed here would matter less, as the strategy could be adjusted to obtain a pattern of aerosol distribution functional to obtaining a certain pattern of cooling. The consistency in the normalized temperature response between models indeed seems to suggest this is possible. If similar large-scale patterns of cooling can be achieved consistently between models, more focus can then be given to understanding the uncertainties in the projected regional-scale responses, which are essential to properly assess local risks and adaptation strategies. This does not remove uncertainties related to the large-scale transport of the aerosols or those arising from any indirect drivers of the surface responses, like the associated changes in concentrations of lower-stratospheric ozone and water vapor alongside any changes in circulation; but it suggests that, with this strategy, they can be more easily separated, as we illustrate in more depth in Part 2 using these same simulations. Due to the potential for our proposed experiments to shed further light on processes relevant for SAI and as a fundamental experiment to develop an intercomparison of more comprehensive SAI strategies, we would like to also suggest them as a possible test-bed experiment for the Geoengineering Model Intercomparison Project.