Upper tropospheric ice sensitivity to sulfate geoengineering

Aside from the direct surface cooling that sulfate geoengineering (SG) would produce, investigations of the possible side effects of this method are still ongoing, such as the exploration of the effect that SG may have on upper tropospheric cirrus cloudiness. The goal of the present study is to better understand the SG thermodynamical effects on the freezing mechanisms leading to ice particle formation. This is undertaken by comparing SG model simulations against a Representative Concentration Pathway 4.5 (RCP4.5) reference case. In the first case, the aerosol-driven surface cooling is included and coupled to the stratospheric warming resulting from the aerosol absorption of terrestrial and solar near-infrared radiation. In a second SG perturbed case, the surface temperatures are kept unchanged with respect to the reference RCP4.5 case. When combined, surface cooling and lower stratospheric warming tend to stabilize the atmosphere, which decreases the turbulence and updraft velocities (−10 % in our modeling study). The net effect is an induced cirrus thinning, which may then produce a significant indirect negative radiative forcing (RF). This RF would go in the same direction as the direct effect of solar radiation scattering by aerosols, and would consequently influence the amount of sulfur needed to counteract the positive RF due to greenhouse gases. In our study, given an 8 Tg-SO₂ yr⁻¹ equatorial injection into the lower stratosphere, an all-sky net tropopause RF of −1.46 W m⁻² is calculated, of which −0.3 W m⁻² (20 %) is from the indirect effect on cirrus thinning (6 % reduction in ice optical depth). When surface cooling is ignored, the ice optical depth reduction is lowered to 3 %, with an all-sky net tropopause RF of −1.4 W m⁻², of which −0.14 W m⁻² (10 %) is from cirrus thinning. Relative to the clear-sky net tropopause RF due to SG aerosols (−2.1 W m⁻²), the cumulative effect of the background clouds and cirrus thinning accounts for +0.6 W m⁻², due to the partial compensation of large positive shortwave (+1.6 W m⁻²) and negative longwave adjustments (−1.0 W m⁻²). When surface cooling is ignored, the net cloud adjustment becomes +0.8 W m⁻², with the shortwave contribution (+1.5 W m⁻²) almost twice as much as that of the longwave (−0.7 W m⁻²). This highlights the importance of including all of the dynamical feedbacks of SG aerosols.

1 Introduction

Sulfate geoengineering (SG) is one of the methods that have been proposed by the scientific community (Budyko, 1977; Crutzen, 2006; Niemeier and Tilmes, 2017) to cool our planet for a limited amount of time, in response to the warming caused by the increasing greenhouse gases (GHGs) of anthropogenic origin. SG proposes the injection of SO₂ into the tropical lower stratosphere in order to produce an optically active cloud of H₂SO₄–H₂O supercooled liquid aerosols that would reflect part of the incoming solar radiation back to space. However, these aerosols would simultaneously warm the lower stratosphere by a few degrees. The idea stems from the cooling effect of past explosive volcanic eruptions in the tropical region (the last being Pinatubo in 1991). These major eruptions injected large amounts of SO₂ into the lower stratosphere and increased the planetary albedo. The resulting cooling effect was clearly observed (Robock, 2000), although its magnitude is still being discussed (Canty et al., 2013).

In the case of past volcanic eruptions, both the direct and indirect effects of episodic large injections of sulfur into the stratosphere have been observed and documented; this is obviously not possible for planned sustained sulfur injections in SG experiments. Because of this, the scientific community mainly relies on simulations using climate models and comparisons of the results among them, such as the GeoMIP project (Kravitz et al., 2011, 2013). Different injection scenarios have been proposed and adopted in modeling experiments, the most common of these being a scenario involving a constant sulfur injection rate at the Equator for a certain number of years to understand the climate response to such an atmospheric perturbation. Simulations have also been performed to identify the optimal magnitude and location of the stratospheric sulfur injection and to obtain the highest ratio between the radiative forcing (RF) and the injection magnitude (Niemeier and Schmidt, 2017; Tilmes et al., 2017; Kleinschmitt et al., 2018).

Amongst the various side effects of SG, those with non-negligible impacts on the RF were analyzed and summarized in Visioni et al. (2017a). These side effects were related to an enhancement of stratospheric ozone destruction (Tilmes et al., 2008; Pitari et al., 2014; Xia et al., 2017), an increase in the concentration and lifetime of methane (Visioni et al., 2017b), an increase of stratospheric water vapor due to a tropical tropopause layer (TTL) warming (Pitari et al., 2014) and, most importantly, to a change in the probability of the formation of cirrus ice particles in the upper troposphere (UT) (Kuebbeler et al., 2012). Regarding this latter effect, some studies have appeared in the recent literature that propose ways in which SG could affect the UT cirrus ice number density and optical depth. We will discuss these studies below and try to expand on certain aspects in the present work.

In an unperturbed atmosphere, the formation of UT ice particles may take place due to either homogeneous or heterogeneous freezing (Karcher and Lohmann, 2002; Hendricks et al., 2011), with the former process normally dominating over the latter, at least in model simulations (Storelvmo and Herger, 2014; Gasparini and Lohmann, 2016; Gasparini et al., 2017; Barahona et al., 2017). However, Cziczo et al. (2013)reported that, in some areas, in situ measurements show that heterogeneous freezing dominates over homogeneous freezing. Homogeneous freezing takes place when the ice saturation ratio is relatively high (typically above ∼1.5), local temperatures are below the threshold for atmospheric ice particle formation (∼238 K) and supercooled solution droplets are present, namely sulfate aerosols or sulfate-coated aerosols. Supersaturation conditions are maintained by intense vertical motions controlling the adiabatic cooling rate and bringing water vapor from the lower to the upper troposphere. Ice crystals formed in this way both reflect part of the incoming solar radiation (negative RF) and trap part of the outgoing planetary radiation, contributing to the greenhouse effect (positive RF). The sign of the combined effects is not easily determined in a variety of atmospheric conditions. Generally, it has been shown that the net UT ice contribution to the RF is positive (Chen et al., 2000; Fusina et al., 2007; Gasparini et al., 2017). However, this is a rather delicate balance and strongly depends on the humidity, cloud cover and optical properties (Mitchell et al., 2008); therefore, a robust atmospheric perturbation, such as the one that the SG could produce, may significantly affect it.

The perturbation to the UT ice could be twofold. On one hand, Cirisan et al. (2013) studied how the H₂SO₄–H₂O droplets resulting from the sulfur injection would interact with cirrus clouds, both microphysically and radiatively. An upper tropospheric increase of the sulfate aerosol number concentration is expected under SG conditions due to gravitational sedimentation and the large-scale transport of particles below the tropopause from the lower stratosphere (LS). However, sulfuric acid liquid supercooled droplets cannot act as ice nuclei (IN) for heterogeneous freezing. At the same time, the background number concentration of the UT aerosols acting as nuclei for homogeneous freezing is already much higher with respect to the ice particle number density. Therefore, a negligible increase of the active IN population would be found in the UT (mainly due to a shift in the distribution of sulfate particles towards radii where homogeneous freezing is more favorable); the same would hold true for the positive RF associated with a possible increase of ice particles from this effect (Cirisan et al., 2013).

Kuebbeler et al. (2012), on the other hand, analyzed the effects produced by dynamical changes due to the modification of the tropospheric thermal gradient produced by stratospheric geoengineering aerosols. In particular, the LS warming, caused by increasing heating rates in the optically thick sulfate cloud, tends to decrease the tropospheric lapse rate. A subsequent decrease in the available turbulent kinetic energy (TKE) would follow and translate in a slowing down of the updraft and the adiabatic cooling rate, which would reduce the probability of sufficiently high supersaturation values capable of producing ice crystals formation via homogeneous freezing. Their study also found a resulting large reduction in the magnitude of the net RF with respect to clear-sky conditions, where only direct aerosol forcing is considered (−0.93 W m⁻² against −1.53 W m⁻²). They concluded that this forcing reduction results not only from the mere (passive) presence of background clouds that affect the atmospheric radiative transfer, but also from the cirrus cloud thinning produced by the SG aerosols. This may have clear implications regarding the potential use of SG to counterbalance global warming.

However, the aforementioned study lacked an important part of the possible dynamical feedback of SG – the changes in sea surface temperatures (SSTs) that would result from the decreased incoming solar radiation. The goal of this research was to study the impact of a stratospheric sulfate injection on cirrus ice particles formed via homogeneous freezing and to understand how both the local stratospheric warming and the surface and tropospheric cooling can affect this process; to accomplish this, we utilized the composition–climate coupled model developed at the University of L'Aquila (ULAQ-CCM). We performed a SG simulation with an 8 Tg-SO₂ yr⁻¹ injection, using surface temperatures (T_s) calculated using the CCSM-CAM4 (Community Climate System Model – Community Atmosphere Model version 4) atmosphere–ocean coupled model, operated with the same sulfur injection (thus resulting in a general surface cooling with respect to unperturbed atmospheric conditions). This perturbed experiment (named G4, according to the convention of Kravitz et al. (2011), regardless of the time constant magnitude of the injection) is compared against a baseline simulation without SG and using a background anthropogenic emission scenario corresponding to the Representative Concentration Pathway 4.5 (RCP4.5) (Taylor et al., 2012) (named “Base” in our study). To properly compare our results with those of Kuebbeler et al. (2012), a third simulation was performed with the same geoengineering sulfur injection as G4 but with the surface temperatures fixed at the Base case values (named G4K).

The effects of the SG T_s changes on the lower stratospheric dynamics have already been discussed in Visioni et al. (2017b); in this study, we focus on their impact in the upper troposphere. Unlike the other side effects of sulfur injection into the stratosphere, a comparison between the effects of a volcanic eruption and SG on cirrus ice is difficult to draw. This is mainly because during a volcanic eruption episode (contrary to SG), a large amount of solid ash particles is injected into the lower stratosphere together with SO₂. Part of these particles, after settling down below the tropopause, may contribute to increasing the number density of IN available for heterogeneous freezing in the UT. This could help to explain some observed increases in UT ice particles after the Pinatubo eruption (Sassen et al., 1995). More recently, Friberg et al. (2005) showed that cirrus cloud reflectance and optical depth are reduced in the Northern Hemisphere in periods with more pronounced volcanic activity. However, other studies such as Meyer et al. (2015) dispute this effect, and no conclusive answer can be given.

Understanding the RF contribution of the UT ice perturbation in a SG scenario is particularly crucial if the scientific community wants to design experiments where the goals are to meet a given climate target, as proposed in Kravitz et al. (2017) and MacMartin et al. (2017).

This paper is structured in three subsequent sections plus the conclusions. In Sect. 2, we describe the CCSM-CAM4 and ULAQ-CCM models and the setup of the numerical experiments; we also try to evaluate the ULAQ-CCM skill in simulating the formation of the cirrus ice clouds, using reanalysis and satellite data. In Sect. 3, we discuss the model-calculated changes in the thermodynamical properties of the atmosphere and in cirrus cloudiness (size distribution, extinction, optical depth, number concentration) produced by the SG. Finally, we also show how these perturbations translate into tropopause radiative forcing terms.

2 Model descriptions and setup of numerical experiments

2.1 CCSM-CAM4

The Community Climate System Model – Community Atmosphere Model version 4 (CCSM-CAM4) is an atmosphere–ocean coupled model that was used in this experiment to calculate the evolution of T_s for both the Base case (RCP4.5 scenario) and a geoengineering case with the same sulfur injection as the ULAQ-CCM model, described in Tilmes et al. (2015). For these simulations, the model was run without interactive chemistry. The resolution of the model is $\mathrm{1.9}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ with 26 vertical levels and the top of the model is at 3 hPa. The model has been fully described in Neale et al. (2013) and Tilmes et al. (2016) and has been shown to compare well against observations in the stratosphere from Lamarque et al. (2012). Ice clouds are diagnosed from a purely relative humidity-based formulation (Neale et al., 2013). The results of an 8 Tg-SO₂ yr⁻¹ injection on surface temperatures and the effects of the inclusion of the perturbed T_s in the ULAQ-CCM model have already been discussed in Visioni et al. (2017b).

2.2 ULAQ-CCM

The first version of the University of L'Aquila composition–climate coupled model was described in Pitari et al. (2002); subsequent model versions have been documented in modeling intercomparison campaigns (Eyring et al., 2006; Morgenstern et al., 2010, 2017). Model updates of the horizontal and vertical resolution, photolysis cross sections, the treatment of Schumann–Runge bands and the radiative transfer code have been described and tested in Pitari et al. (2014) and Chipperfield et al. (2014). The shortwave radiative module has been documented and tested for tropospheric aerosols in Randles et al. (2013) and for volcanic stratospheric aerosols in Pitari et al. (2016a). It makes use of a two-stream delta-Eddington approximation and is online in the model for photolysis, solar heating rates and radiative flux calculations. A companion broadband, k distribution longwave radiative module is used for the heating rate and radiative flux calculations in the planetary infrared spectrum (Chou et al., 2001).

A critical atmospheric region in SG studies is the upper troposphere–lower stratosphere (UTLS). An extensive model evaluation based on specific physical and chemical aspects was made in Gettelman et al. (2010) and Hegglin et al. (2010). Subsequent model improvements in this region were discussed in Pitari et al. (2016b). The treatment of T_s, and their importance for the lower stratospheric dynamics and species transport under a geoengineering scenario, has been discussed in Visioni et al. (2017b). Another very important aspect to be taken into account for large-scale species transport in the lower stratosphere is the role of the quasi-biennial oscillation (QBO) in SG studies. This has been discussed from different points of view in some recent studies (Aquila et al., 2014; Niemeier and Schmidt, 2017; Visioni et al., 2018). A nudging procedure for the QBO is adopted in the ULAQ-CCM, based on an observed historical data series of equatorial mean zonal winds (Morgenstern et al., 2017).

Table 1 Summary of ULAQ-CCM features and numerical experiments for the present study.

^a Sea surface temperatures from observations; online explicitly calculated land temperatures. ^b Surface temperatures from CCSM-CAM4 (land, ocean, sea ice coverage), separately for Base and G4 (Visioni et al., 2017b); Base values also used for G4K. Indirect effects of SG aerosols on surface temperatures are calculated online in the ULAQ-CCM radiative module (due to UT ice, GHGs and SO₄ imbalance relative to CCSM-CAM4); see text in Sect. 2.3. “eqt. wind obs.” refers to equatorial wind observations.

For the sake of completeness, we discuss some of the model features (fully illustrated in Table 1) in the following, in particular, those relevant for stratospheric sulfate aerosols and upper tropospheric cirrus ice particle formation.

2.2.1 Stratospheric sulfate aerosols

In SG experiments G4 and G4K, SO₂ is injected at the Equator (0^∘ longitude) throughout the altitude range from 18 to 25 km with a Gaussian distribution centered at 21.5 km. The OH oxidation of SO₂ starts the production of supercooled H₂O-H₂SO₄ particles; the size distribution of these particles is calculated in an aerosol microphysics module with a sectional approach, starting from gas–particle interaction processes (nucleation, H₂SO₄ condensation and H₂O growth) and then including aerosol particle coagulation. Removal processes are included via gravitational settling across the tropopause and evaporation in the upper stratosphere (Visioni et al., 2018).

In the troposphere, the ULAQ-CCM model includes sulfate production from dimethyl sulfide (DMS) and SO₂ emissions, with gas phase and aqueous/ice SO₂ oxidation (by OH and H₂O₂, O₃, respectively) to produce SO₄ (Feichter et al., 1996; Clegg and Abbatt, 2001). The tropospheric and stratospheric SO_x budget in the ULAQ-CCM (for unperturbed background conditions) was recently discussed in Pitari et al. (2016c), with a focus on the role of non-explosive volcanic sulfur emissions, and in Visioni et al. (2018), in connection with the SG.

Aerosol extinction, optical thickness, single scattering albedo and surface area density are calculated online at all model grid points every hour. This allows the interactive calculation of up/down diffuse radiation and absorption of solar near-infrared and planetary radiation by SG aerosols, with explicit full coupling of the aerosol, chemistry and radiation modules in the ULAQ-CCM model. This justifies the “composition–climate” name for this coupled model, which is more general than the usual “chemistry-climate” model denomination.

The ability of ULAQ-CCM to produce the correct confinement of sulfate aerosols in the tropical stratosphere has already been documented in the literature. This was carried out by a comparison with SAGE II data following the Pinatubo eruption or looking at the SG conditions during the period following the eruption (see Pitari et al., 2014, 2016a; Visioni et al., 2017b).

2.2.2 Upper tropospheric ice

The formation of UT ice particles may take place via heterogeneous and homogeneous freezing mechanisms. In the latter case, the ULAQ-CCM model adopts the approach initially described in Karcher and Lohmann (2002), which assumes ice crystals are only formed via the homogeneous freezing of solution droplets as a function of local UT temperatures and updraft velocities; this approach also includes the effects of a variable aerosol size distribution. These updraft velocities are obtained as the sum of a dominant term related to the TKE and a much smaller contribution from the large-scale tropospheric circulation (Lohmann and Karcher, 2002). Typical vertical velocity net values are on the order of 10–20 cm s⁻¹ (see Sect. 3.1) and allow the formation of thin cirrus.

For the ice supersaturation ratio, we adopt a simplified probabilistic approach, starting from the knowledge of climatological frequencies of the UT relative humidity (RH_ICE), from which a mean value and a standard deviation can be calculated, assuming a normal distribution. Local ice supersaturation conditions (RH_ICE>100 %) are a result of turbulent ascent and can be found in the UT, in the vertical layer below the tropopause (where turbulent updraft conditions may be found) and above an altitude where T<238 K (i.e., the assumed threshold for the spontaneous freezing of solution droplets). Here, the conditions for ice formation are met and we may calculate the probability that RH_ICE>1.5 (P_HOM). This represents the assumed threshold for the activation of homogeneous freezing (in our model this threshold does not depend on local temperature or water activity conditions), which is considerably higher with respect to the threshold for heterogeneous freezing (RH${}_{\mathrm{ICE}}>\sim \mathrm{1.3}$) (Hendricks et al., 2011). This represents the probability that an ice particle could be formed via heterogeneous freezing (P_HET) on a preexisting population of ice condensation nuclei (N_IN), typically mineral dust or black carbon (BC) particles transported from the surface.

The size distribution and number density of ice particles formed via heterogeneous freezing (n_HET) is calculated utilizing the formulation of Hendricks et al. (2011) and the ULAQ microphysical scheme adopted for polar stratospheric ice particle formation (Pitari et al., 2002). N_IN is the sum of grid point model-predicted concentrations of mineral dust (DU) and BC aerosols (N_DU and N_BC, respectively) and is used as the population of available condensation nuclei, with P_HET being the probability that RH_ICE>1.3 at any model grid point. The problem in this case is the actual availability of solid ice nuclei. A low fraction of activated IN is suggested in the literature (f_DU=1 % for mineral dust and f_BC=0.25 % for BC) because the large majority of IN will be rapidly coated by sulfate (Hendricks et al., 2011). The number density n_HET is then obtained as

$$\begin{array}{}\text{(1)}& n_{\mathrm{HET}}=\left(f_{\mathrm{BC}}{N}_{\mathrm{BC}}+f_{\mathrm{DU}}{N}_{\mathrm{DU}}\right)P_{\mathrm{HET}}.\end{array}$$

The specification of the active ice fraction for both mineral dust and BC represents the major source of uncertainty for UT ice particle formation via heterogeneous freezing. Considering the above assumptions, homogeneous freezing normally dominates ice particle formation in the ULAQ-CCM model, with respect to the heterogeneous freezing mechanism. Whilst this may not be considered a general conclusion, it is assumed to be valid in all thermodynamics conditions and any local atmospheric composition; this has been shown, for instance, in Cziczo et al. (2013), where a predominance of heterogeneous freezing over homogeneous was found. In general, the freezing mechanism that dominates in the atmosphere is still very uncertain.

The calculated mass mixing ratio of ice formed in the ULAQ-CCM model through both freezing mechanisms is shown in Fig. 1a, c for two pressure layers, 150–200 and 350–400 hPa, where the ice formation is greater in the tropics and mid–high latitudes, respectively. These calculations are compared against the MERRA-2 (Bosilovich et al., 2017; Gelaro et al., 2017) and ERA5 reanalyses (Stephens et al., 2000), and are all averaged over the same decade (2003–2012). For the upper layer (150–200 hPa), we also show the MLS satellite retrieval in Fig. S1 in the Supplement (Wu et al., 2008), which compares very closely to the ERA5 reanalysis. Tropical ice formation shows a strong land–ocean asymmetry due to significantly higher P_HOM and P_HET values over land. For both pressure layers, the magnitude and spatial distribution of the ice mass mixing ratio are comparable between the ULAQ-CCM and MERRA. Regarding the datasets used to compare against our model results, note that there is a large spread amongst retrievals (such as MODIS or CALIPSO) and amongst reanalyses (Zhang et al., 2010; Duncan and Eriksson, 2018). In particular, MERRA-2 appears to be at the lower end of the spectrum with regards to some quantities, such as ice water path. Considering that the dataset only considers non-precipitating ice (Duncan and Eriksson, 2018), this quantity might be closer to the one simulated in our model and thus would allow for a more correct comparison.

Figure 1 Lat–long maps of the ice mass mixing ratio (mg kg⁻¹-air) for pressure layers representative of the tropical (a, c, e) and extratropical (b, d, f) upper troposphere. Panels (a) and (b) are for ULAQ-CCM; panels (c) and (d) are for MERRA-2 data (Bosilovich et al., 2017); panels (e) and (f) are for ERA5 data (Stephens et al., 2000). Time is averaged over the years 2003–2012.

While the probability of homogeneous ice formation is defined as above, the number density and size of the ice particles formed this way is determined by the local temperatures and vertical velocities, in addition to the competing ice formation mechanism – heterogeneous freezing. The lower the temperature, the faster the nucleation rate; thus, more ice crystals can be formed. Conversely, higher vertical velocities increase the saturation ratio, leading to more ice crystals being formed before water deposition on ice crystals reduces supersaturation below the threshold. The spatial distribution of the cirrus ice optical depth (OD) in the model is calculated as follows:

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{\mathrm{ice}}=\mathrm{\Delta }z\sum _i\sum _j{Q}_{\mathrm{ext}}\mathit{\pi }{r}_{ij}^{\mathrm{2}}{n}_{ij}\left(r\right),\end{array}$$

where the extinction efficiency coefficient Q_ext∼2 at all visible wavelengths for ice particle sizes is on the order of 5–50 µm; i is an index for the vertical layers, and the sum is over all the vertical layers in the UT; j is an index for the particle size bins, and the sum is over the whole size distribution; r_ij is the particle radius at the ith layer and jth bin; and n_ij is the corresponding ice number density.

Equation (2) can easily be applied to the model, and the results are shown in Fig. 2a. An evaluation can again be made using the ice mixing ratio from MERRA-2 and ERA5 (shown in Fig. 1b, c for two specific pressure layers), and the ULAQ-CCM values of the ice particle effective radius. With these two quantities we indirectly derived τ_ice at every horizontal grid point in Eq. (2), using the hydrostatic equation:

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{ice}}=Q_{\mathrm{ext}}{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{\Delta }p}{g}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{ice}}}}\sum _i{\displaystyle \frac{{\mathit{\chi }}_i}{r_i}},\end{array}$$

where the sum is, again, over all the vertical layers (constant Δp=50 hPa), g is the acceleration of gravity, ρ_ice is the ice bulk density, r_i is the ULAQ-CCM effective radius at the ith layer, and χ_i is the MERRA-2 and ERA5 ice mass mixing ratio at the ith layer. Through this process, we obtain the optical depth values in Fig. 2b, c. The ODs are comparable in terms of spatial distribution, with the highest values in the tropics over land. However, the absolute values in the ULAQ-CCM model are significantly smaller over the tropics. The reason for this is that updraft velocities result in a relatively narrow interval (w<30 cm s⁻¹) when only calculated as a function of TKE (as in the ULAQ-CCM), while thick cirrus formation takes place from strong (and less frequent) convective events (w<100 cm s⁻¹). This detrained ice that originates in deep convection is not included in our model formulation.

Figure 2 Ice optical depth at λ=0.55 µm, from ULAQ-CCM calculations (a) and from the MERRA-2 (b) and ERA5 (c) ice mass mixing ratio ($\mathrm{100}<p<\mathrm{450}$ hPa), with ULAQ-CCM particle effective radius. Time is averaged over the years 2003–2012.

Figure 3 Fraction of total ice formed through heterogeneous freezing in ULAQ-CCM averaged over the years 2003–2012, as a function of latitude and longitude for the total optical depth (a) and as a function of altitude and latitude for the zonally averaged extinction (b). In panel (b), the color scale is logarithmic, ranging from 0.01 (i.e., 1 % of total ice extinction) up to 1 (100 %). The dashed lines show the mean tropopause height with seasonal variability (where seasonal variability is defined as ±1σ of the average height). The dash-dotted lines show the mean height (with seasonal variability) at which T=238 K (freezing is allowed for colder temperatures).

In Fig. 3, we show the model-predicted fraction of ice formed through heterogeneous freezing in terms of optical depth (Fig. 3a) and zonally averaged extinction (Fig. 3b). In both panels, we see that a large part of the ice particles formed via heterogeneous freezing is located in the tropical band at lower altitudes, where a higher concentration of mineral dust and BC ice nuclei can be transported from the surface. In these regions, the fraction of ice formed this way can be as much as 80 % of the total.

In Fig. 4a, b we show the model calculated vertical profiles of the ice particle number density averaged over the tropics (Fig. 4a) and the extratropics (Fig. 4b), superimposed with the time variability produced by changing conditions of vertical velocity, temperature, and P_HOM and P_HET. The ice number density maxima are located at different altitudes in the two latitude bands, close to 13 km in the tropics and 8 km elsewhere. This is clearly expected from the latitudinal variability of the tropopause height.

With a procedure similar to the one described above for the ice OD, we may derive a first-order approximation of the ice number density from the MERRA-2 and ERA5 ice mass mixing ratio and ULAQ-CCM radii. Similar to Eq. (3), for the ice number density (n_i) at each vertical layer we obtain the following expression:

$$\begin{array}{}\text{(4)}& n_{\mathrm{i}}={\displaystyle \frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{ice}}}}{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{3}}}}{\mathit{\chi }}_i.\end{array}$$

Figure 4 Average upper tropospheric profiles of the ice particle number density (cm⁻³), for the tropics (25^∘ S–25^∘ N; panel a) and extratropics (35–90^∘ S, 35–90^∘ N; panel b). Time is averaged over the years 2003–2012. The shaded areas represent ±1σ for the ensemble over the 10-year period. The red circles show indirectly derived values from the MERRA-2 and ERA5 ice mass mixing ratio and ULAQ-CCM effective radius (see text). Panel (c) shows the zonally averaged probability of ice formation via homogeneous freezing (percent), as a function of altitude and latitude. The dashed lines show the mean tropopause height (with seasonal variability). The dash-dotted lines show the mean height at which T=238 K (with seasonal variability).

The results from Eq. (4) (circles in Fig. 4a, b) show that while the model and the indirectly derived values from the reanalyses agree in terms of the general vertical distribution and localization of the maxima in the extratropics, the ULAQ-CCM tends to have smaller number densities in the tropics in the 10–13 km layer. Again, this should not be surprising in light of the fact that we are focusing on a specific type of cirrus cloud particle.

Table 2 Summary of globally and time-averaged sulfate aerosol and cirrus ice particle related quantities, as calculated in the ULAQ-CCM model and compared with available satellite and reanalyses data. Sulfate aerosols – sectional approach (Pitari et al., 2002, 2014). Cirrus ice particles: parameterization for homogenous (HOM) and heterogeneous (HET) freezing are summarized in the text and based on the formulation from Karcher and Lohmann (2002) (HOM), but including the effects of the aerosol size distribution, and Hendricks et al. (2011) (HET); a probabilistic approach is adopted for the ice supersaturation ratio. Standard deviations are calculated over the time series of globally averaged monthly mean values. On the global average, our model predicts a 90 % fraction of the ice optical depth formed via homogeneous freezing.

Figure 4c shows the model-calculated values of P_HOM, as a 2-D zonally averaged distribution. Using these P_HOM values, it is possible to scale a n_i value measured in the midlatitude airborne campaign of Strom et al. (1997) during a young cirrus formation, to derive an average climatological value to be considered consistent with our modeling approach. Strom et al. (1997) measured a midlatitude ice concentration value of n=0.3 cm⁻³ in a young cirrus cloud at T=220 K and p=320 hPa. If we scale this result with our corresponding $P_{\mathrm{HOM}}=\mathrm{12}\pm \mathrm{3}$ %, a “climatological-mean” value $n=\mathrm{0.025}\pm \mathrm{0.005}$ cm⁻³ is obtained, which is close to our model prediction value of 0.031±0.008 cm⁻³ (Fig. 4b).

Relevant aerosol and ice quantities calculated in the ULAQ-CCM model are summarized in Table 2 and compared with available satellite observations. The first two rows in Table 2 compare the ULAQ-CCM results for stratospheric sulfate optical depth (OD) and the tropical effective radius (r_eff) against SAGE-II and AVHRR satellite observations (Thomason et al., 1997; Long and Stowe, 1994), under post-Pinatubo conditions (Pitari et al., 2016a). This is undertaken to highlight the realistic representation of the gas–particle conversion and aerosol microphysics processes in the model, along with the aerosol large-scale transport in the lower stratosphere in the case of a major tropical volcanic eruption, which may be used as a proxy for SG with an equatorial SO₂ injection. A comparison of the aerosol effective radii under volcanic and background conditions (see rows 2 and 3 in Table 2) clearly shows the effects of the sulfuric acid condensation on the size extension of the aerosol accumulation mode and how this is represented in the model.

The bottom 5 rows in Table 2 compare the global budget calculations for tropospheric ice particles with values obtained from the MERRA-2 and ERA5 reanalyses (ice mass mixing ratio) and ULAQ-CCM effective radius (compared in row 7 with the ice effective radius as retrieved by MODIS). The simultaneous use of these two products (reanalysis values for ice mass mixing ratio and model calculated radius) allows for an indirect calculation of the ice optical depth (row 8 of Table 2), as previously discussed. The ULAQ-CCM OD underestimation is mostly related to the lower values of the ice mass mixing ratio in the largest portion of the upper troposphere (see row 5 of Table 2) and may be, in part, explained with the inclusion of a relatively narrow interval for updraft velocities (w<30 cm s⁻¹).

The values are given separately for the ice formed through homogeneous and heterogeneous freezing.

2.3 Setup of the numerical experiments and role of perturbed SSTs

The use of a composition–climate coupled model, such as the ULAQ-CCM model, offers multiple advantages in this type of study: (a) the online inclusion of interaction between aerosol and ice particles microphysics with chemistry, radiation, climate, dynamics and transport; (b) the stratosphere–troposphere explicit interactions for the large-scale transport of gas and aerosol species (the model adopted high vertical resolution is important across the tropopause region); (c) the sufficiently detailed chemistry both the stratosphere and troposphere, with a robust design for heterogeneous chemical reactions on sulfuric acid aerosols, polar stratospheric cloud particles, and upper tropospheric ice and liquid water cloud particles. This allows us to account for the atmospheric circulation changes produced by sulfate geoengineering. The ULAQ-CCM model has proven many times to be capable of producing sound physical and chemical responses to both sulfate geoengineering (Pitari et al., 2014; Visioni et al., 2017b) and for large explosive volcanic eruptions (Pitari, 1993; Pitari et al., 2016a, b).

Figure 5 Seasonally averaged surface temperature anomalies G4-RCP4.5 (K) from the CCSM-CAM4 atmosphere–ocean coupled model (time-averaged over the years 2030–2069). The shaded areas are not statistically significant within ±1σ. Panels (a)–(d) refer to the following periods: December–January–February (a); March–April–May (b); June–July–August (c); and September–October–November (d).

In addition to a reference historical model experiment (1960–2015), we performed three sets of SG simulations: a baseline (Base) unperturbed case and two geoengineering experiments (G4 and G4K), both run with an injection of 8 Tg-SO₂ yr⁻¹ into the equatorial stratosphere between 18 and 25 km of altitude. A very similar experimental design was described in Kravitz et al. (2011) for the GeoMIP G4 experiment with a sustained fixed injection of sulfur dioxide (5 Tg-SO₂ yr⁻¹ in that case; we use 8 Tg-SO₂ yr⁻¹, although all other prescriptions such as height and the latitude of the injections are the same as in the abovementioned paper). These numerical experiments were all run between the years 2020 and 2069, with analyses focusing on the 2030 to 2069 period; all experiments also take place under the same RCP4.5 reference scenario for well-mixed greenhouse gases. The ULAQ-CCM model is not an atmosphere–ocean coupled model and uses externally provided surface temperatures as prescribed boundary conditions for the dynamical module. These surface temperatures are taken from the CCSM-CAM4 model, which was run under the same RCP4.5 and G4 conditions (8 Tg-SO₂ yr⁻¹ fixed injection into the equatorial lower stratosphere). In this way our main experiment G4 may account for the T_s response to SG (Fig. 5). We acknowledge that this procedure may only be valid as a first-order approximation, considering that CCSM-CAM4 was not run with a coupled chemistry and utilized a much simpler cirrus parameterization that produces negligible changes in the geoengineering experiment (Neale et al., 2013). However, we believe it to still be a consistent one, considering that the main effect produced by the sulfate injection is the direct aerosol effect (Visioni et al., 2017a), and that the prescribed stratospheric aerosol field in the SG simulation in CCSM-CAM4 (Tilmes et al., 2015) is comparable to the one produced by the sulfate injection in ULAQ-CCM. With this in mind, in the next paragraphs, we first discuss the T_s perturbation and its significance for this study and then the approach adopted for minimizing the inconsistency introduced in ULAQ-CCM due to the use of T_s from a different model.

A strong inter-hemispheric asymmetry in the T_s changes produced by SG is evident in Fig. 5 (see also the annually and zonally averaged values in Fig. 6a), with a negative anomaly in the Arctic region that is approximately 1 K larger than that of the high southern latitudes. The SG cooling impact on the Arctic sea ice is such that larger negative surface temperature anomalies are favored in the Northern Hemisphere high latitudes for several months during the year, from the fall to spring months (see Fig. 5a, b, d); thus, atmospheric stabilization is increased in the Northern Hemisphere with respect to the Southern Hemisphere. Note, however, that the dynamical effects of this enhanced atmospheric stability in the SG conditions (decreasing wave activity and turbulence) may be partially counterbalanced by the increased longitudinal variability of the induced cooling, mostly connected with positive surface temperature anomalies in the subpolar North Atlantic. These positive temperature anomalies in the North Atlantic sub-Arctic are a direct consequence of the increasing amount of polar sea ice in the SG conditions, with the southward transport of colder and saltier ocean waters in the sub-Arctic, compared with RCP4.5 Base conditions (Tilmes et al., 2009). In this way, the North Atlantic subpolar downwelling of these cold surface waters to the deep ocean is favored with respect to the Base conditions, which produces positive anomalies in sea surface temperatures.

Figure 6 (a) Zonally averaged T_s anomalies G4-RCP4.5 (K), under different conditions for the G4 perturbed case (time-averaged over the years 2030–2069): from the atmosphere–ocean coupled model CCSM-CAM4 (black dashed line); as above, but adding the T_s anomalies in the ULAQ-CCM model with online coupling of cirrus ice changes and the SO₄ imbalance between CCSM-CAM4 and ULAQ-CCM (green line); as above, but also adding the T_s anomalies in the ULAQ-CCM model with online coupling of GHG changes (red line) (see text and legend). The shaded area represents ±1σ of the zonally averaged T_s over the 40-year period. (b) Lat–long distribution of the T_s anomalies (K) calculated online in the ULAQ-CCM model considering cirrus ice changes, the SO₄ imbalance between CCSM-CAM4 and ULAQ-CCM and GHG changes (time-averaged over the years 2030–2069).

Although not statistically significant, the SG-induced warming on the Antarctic continent during wintertime (Fig. 5c) is a direct consequence of the geoengineering aerosol positive radiative forcing in the planetary longwave, which represents the net forcing at these high latitudes in the absence of sunlight. This radiative feature will be further discussed in Sect. 3. All of these high-latitude positive temperature anomalies directly reflect in the large variability of the zonally averaged surface temperature changes presented in Fig. 6a.

To correct for the potentially significant model inconsistency introduced due to the use of surface temperatures taken from a different model, the following procedure has been adopted. The ULAQ-CCM radiative-climate module has been modified for calculating the T_s perturbation produced by the radiative flux changes due to the stratospheric sulfate aerosol imbalance online (in a fully coupled approach) with respect to the CCSM-CAM4 distribution in the G4 case. In addition, we also include the SG-driven indirect perturbation of greenhouse gases (O₃, H₂O, CH₄ and CO₂ from the changing methane oxidation) in the radiative balance, as well as the indirect perturbation of upper tropospheric ice particles. This online calculated T_s perturbation is then added to the externally provided T_s field from CCSM-CAM4 for the G4 experiment. Table S1, Figs. S2 and 6 document these radiative flux changes and their impact on the calculated T_s.

Surface temperature changes due to the abovementioned indirect SG effects are calculated from the instantaneous perturbation of radiative fluxes, which is an exact procedure over continents and polar ice caps, whereas it is only approximate over the oceans. In contrast, as explained above and clearly visible in Table S1, Figs. S1 and 6, the addition from the radiative perturbation to the dominant perturbation (i.e., the one produced by stratospheric sulfate aerosols in the CCSM-CAM4 simulation) is normally small, both globally and locally (notice the different color scale between Figs. 6b and 5). Only the ice induced changes of T_s may be comparable in magnitude to those from the stratospheric aerosols, but are limited to tropical continental surfaces, where UT ice may have significant optical depth values. Furthermore, the SST calculated changes due to chemistry and ice indirect effects of SG are usually smaller, meaning that the impact of our approximation may be negligible.

A sensitivity case (G4K) was run in combination with the G4 simulation, with surface temperatures fixed at the RCP4.5 Base values. Here, the experimental approach is similar to that of Kuebbeler et al. (2012) who ran a G4 simulation with a 5 Tg-SO₂ yr⁻¹ injection and prescribed sea surface temperatures and sea ice from the RCP4.5 Base case. This is undertaken not only to highlight the role of the tropospheric temperature perturbations in cirrus ice formation (given a certain vertical velocity change), but mostly to calculate the updraft sensitivity to different tropospheric stabilization conditions introduced by the stratospheric sulfate aerosol injection.

3 Model response to sulfate geoengineering

In this section, we will show the ULAQ-CCM response to the stratospheric sulfate injection. Some of the perturbations have already been discussed in previous works, in particular those regarding stratospheric dynamics changes (Pitari et al., 2014; Visioni et al., 2017b). Here, we will focus on the thermodynamical changes in the upper troposphere and, consequently, on changes in the formation of cirrus ice clouds.

3.1 Thermodynamical changes in the troposphere

Figure 7 shows the differences in temperature and updraft in G4 and G4K with respect to the Base case. In G4, we observe a tropospheric cooling of ≃1–2 K in the ice formation region throughout all latitudes, while the warming due to the sulfate aerosol absorption of shortwave and longwave radiation is confined above the tropopause (Fig. 7a). When surface temperatures are kept fixed at the RCP4.5 baseline values with the SG perturbation (G4K case), the upper troposphere–lower stratosphere temperature anomalies look very different (Fig. 7b). The tropospheric cooling is absent and the stratospheric warming produced by absorption of longwave planetary and near-infrared solar radiation is more uniformly spread across the lower stratosphere, with some penetration also occurring in the UT (≃0–1 K). The latter is owing to the sulfate aerosol cross-tropopause fluxes that are due to the large-scale transport (at midlatitudes) and gravitational sedimentation (mostly relevant in the tropical region).

The updrafts responsible for the upper tropospheric ice particle formation result from the sum of a rather small large-scale vertical velocity contribution (on the order of 1–2 cm s⁻¹) and a dominant contribution due to motions associated with synoptic scale disturbances and gravity waves (on the order of 10–20 cm s⁻¹); the latter is calculated as a function of the TKE (Lohmann and Karcher, 2002) with the exact formulation reported in Eq. 5:

$$\begin{array}{}\text{(5)}& w_{\mathrm{TOT}}=w_{\mathrm{LS}}+\mathrm{0.7}\sqrt{\mathrm{TKE}}.\end{array}$$

The vertical velocity is reduced in G4 with respect to the Base case by ≃1–2 cm s⁻¹ in the whole UT (Fig. 7c) (on the order of −10 %, as visible in Fig. 8), due to the atmospheric stabilization caused by a reduction in the temperature vertical gradient.

Figure 7 Zonally and time-averaged changes of temperature (a, b) and vertical velocity (c, d) in the G4 (a, c) and G4K (b, d) experiments with respect to the Base case (years 2030–2069). The dashed lines show the mean tropopause height (with seasonal variability). The dash-dotted lines show the mean height at which T=238 K (with seasonal variability). The dotted white lines in panel (a) highlight where ΔT=0 K.

Figure 9a shows the average tropical vertical profiles of the SO₄ mixing ratio (in the particulate phase), for both the Base and SG experiments (with an 8 Tg-SO₂ injection). The changes in the zonally averaged net heating rates, temperatures and zonal winds are also shown in Fig. 9b, c and d, respectively. These values help explain how the SG sulfate perturbation may act as driver for dynamical changes in the UT, with significant effects on ice particle formation.

Figure 8 Average upper tropospheric profiles of the vertical velocity (cm s⁻¹) in the G4 and Base experiments (years 2030–2069). Panels (a) and (b) are for the tropics and extratropics, respectively (see legends). The vertical velocity w is obtained as the sum of the large-scale value and that calculated as a function of the TKE (see Lohmann and Karcher, 2002 and Eq. 5), which essentially accounts for the synoptic scale and gravity wave motions. The shaded areas of the same color represent ±1σ for the ensemble over the 40-year period from 2030 to 2069.

In Fig. 9a, it is interesting to note a somewhat smaller tropical aerosol confinement in the G4K case. This is consistent with the findings of Visioni et al. (2017b): the aerosol-driven surface cooling in G4 (contrary to G4K) favors a decreased wave activity and a consequent decrease in poleward mass fluxes from the tropical reservoir, for both gas and aerosol species. Conversely, the increased H₂SO₄ tropical amount available for aerosol formation tends to produce larger particles with smaller equivalent optical thickness (see Niemeier and Schmidt, 2017; Visioni et al., 2018). In light of this, smaller stratospheric heating rate anomalies are calculated in G4 than in G4K (Fig. 9b): in the latter case, we then expect an enhanced temperature increase in the tropical lower stratosphere (Fig. 9c), coupled to a slight tropospheric warming due to the SG aerosol sedimentation below the tropopause. However, the latter is greatly outweighed by mid–upper tropospheric cooling in G4, due to less intense latent heat exchange resulting from the aerosol-driven T_s decrease (contrary to G4K). As a result, the G4 atmosphere is more efficiently stabilized with respect to G4K, and the positive/negative anomalies of the T∕u (temperature ∕ mean zonal winds) shears in the UT (Fig. 9c, d) favor a decrease of the TKE (and updraft velocities) in G4 compared with G4K (Fig. 7c, d).

Figure 9 Average tropical vertical profiles (25^∘ S–25^∘ N; years 2030–2069) of the SO₄ volume mixing ratio for the G4, G4K and Base experiments (ppbv, panel a); G4–Base changes of net, shortwave and longwave heating rates (K day⁻¹, panel b) (LW is calculated with temperature fixed at Base values) (net heating rate changes are also shown for G4K–Base, with the blue line); G4–Base and G4K–Base temperature changes (K, panel c); G4–Base and G4K–Base changes of mean zonal winds (m s⁻¹, panel d). The shaded areas of the same color represent ±1σ for the ensemble over the 40-year period from 2030 to 2069.

All features of the SW and LW heating rate anomalies in Fig. 9b can be fully explained by taking the aerosol–O₃ coupled effects (Pitari et al., 2014) into account. The sign of tropical ozone changes under SG conditions depends on altitude. The O₃ decreases below ∼25 km and increases above this height; this helps explain the positive/negative heating anomalies in the SW and LW components above an altitude of 25 km.

The SG induced reduction of updraft velocities is significantly smaller in the G4K case (≃0.5 cm s⁻¹, on the order of −3 % the baseline values), as clearly visible in Fig. 7d. This represents the major change in our approach to studying the UT ice sensitivity to SG with respect to the approach adopted in Kuebbeler et al. (2012). According to our calculations, when taking both the main radiative effects of geoengineering stratospheric aerosols (i.e., lower stratospheric heating and surface and tropospheric cooling) into account, the resulting impact on tropospheric turbulence and updraft is significantly enhanced with respect to the case in which only the stratospheric warming is considered. A noticeable difference in the G4K w anomalies with respect to those of G4 is present at low altitudes over the polar regions, where the G4K negative values are larger than in G4. This may be primarily explained by the increasing longitudinal variability of surface temperatures in the G4 case, mainly in the sub-Arctic region (see previous discussion relative to Fig. 5).

The tropical and extratropical average profiles of the updraft velocity are shown in Fig. 8 for both the Base and G4 conditions. The G4K curve (not shown) is an intermediate between the previous two. The pronounced variability of the vertical velocity is expected as a consequence of time, latitude and longitude fluctuations of the TKE. This will produce a significant dispersion of the ice particle size distribution (see ahead in Sect. 3.2).

3.2 Tropospheric ice perturbations due to sulfate geoengineering

In Sect. 2.2.2, we showed that the ULAQ-CCM parameterization for ice particle formation through both homogeneous and heterogeneous freezing produces a spatial distribution of the UT ice particles reasonably comparable to available data in terms of ice number concentration, OD, mass mixing ratio and effective radius. We now move to analyze the model-calculated SG perturbation of some of these quantities by comparing the G4 and G4K simulations against the Base case. As we have previously discussed and shown in Figs. 7–9, these perturbations are essentially produced and regulated by decreasing vertical velocities (−1.7 and −0.8 cm s⁻¹, in the tropical region below the tropopause for G4 and G4K, respectively) and by changing the tropospheric temperatures (−1.2 and +0.5 K, in the tropical UT region for G4 and G4K, respectively).

The model-calculated globally and time-averaged size distribution of the ice particles is presented in Fig. 10 for the three experiments, along with their globally averaged effective radius. A significant change in size distribution is highlighted in Fig. 10 in both SG experiments with respect to not only the Base case, but also G4 and G4K. The common feature in both SG cases is the expected decreased particle population over the whole radial spectrum with respect to the Base experiment. This is due to the increased atmospheric stabilization forced by the SG aerosols with reduced updraft velocities and the consequent decrease of the UT ice supersaturation probability.

Figure 10 Globally and time-averaged number density values of ice crystals as a function of particle radius (dn∕dlogr, cm⁻³) (years 2030–2069). Shaded areas of the same color represent ±1σ for the ensemble over the 40-year period from 2030 to 2069. The calculated global mean values of the ice particle effective radius are as follows: Base is  31.2±3.2 µm; G4 is  33.5±3.6 µm; and G4K is  36.8±4.1 µm. The reference MODIS value in Table 2 is 33.4±2.1 µm.

However, the UT temperature anomalies are very different in the two SG experiments with respect to the Base case (see Fig. 7). As a consequence of this, the tropospheric cooling produced in G4 by the T_s adjustment to the stratospheric aerosol negative RF favors a number density increase of ice particles with respect to the G4K experiment but is still less than in the Base case (see also Fig. 11); this is due to the dominant impact of the reduced updraft. Cooler temperatures actually cause a faster nucleation of the ice particles, quickly removing water vapor available for the freezing itself and limiting the condensational growth of ice particles (Kuebbeler et al., 2012; Visioni et al., 2017a). At the same time, the velocity and temperature negative anomalies partially compensate each other also in the particle size spectrum, with a resulting effective radius in G4 larger with respect to the one in the unperturbed atmosphere (33.5±3.2 and 31.2±3.2 µm, respectively) but smaller than that in G4K. In this latter case, the UT is slightly warmed up with respect to the Base case (see Fig. 7) so that both the velocity and temperature anomalies tend to increase the particle size (36.8±4.1 µm). Globally, the ULAQ-CCM baseline values of the effective radius fall well inside the MODIS range of variability (33.4±2.1 µm).

As visible in Fig. 11 the calculated ice number densities follow the zonal mean behavior of the MERRA-2 and ERA5 indirectly derived values, with the previously discussed underestimation tendency, mainly in the tropical region (see Fig. 4).

Figure 11 Zonally and time-averaged total number density values of ice crystals as a function of latitude (n, cm⁻³) (years 2030–2069), as calculated in the ULAQ-CCM (for Base, G4, G4K experiments) and compared with indirectly derived values from the MERRA-2 and ERA5 ice mass mixing ratio and ULAQ-CCM (Eq. 4). Number densities are calculated at pressure layers 150–200 hPa for 25^∘ S–25^∘ N, 200–250 hPa for 25–35^∘ N/S, 250–300 hPa for 35–45^∘ N/S, 300–350 hPa for 45–55^∘ N/S and 350–400 hPa for 55–90^∘ N/S.

Figure 12 Average upper tropospheric profiles of ice particle extinction (λ=0.55 µm) (km⁻¹) for the tropics (25^∘ S–25^∘ N) and extratropics (35–90^∘ S, 35–90^∘ N) in panels (a) and (c) and (b) and (d), respectively. (a, b) Ice extinction changes for G4-Base (red curves) and G4K-Base (blue curves) (years 2030–2069). (c, d) Comparison of ULAQ-CCM calculated values of ice extinction with indirectly derived values from the MERRA-2 and ERA5 ice mass mixing ratio and ULAQ-CCM effective radius (red and blue circles) (see text). The time average is over the years 2003–2012. The shaded areas represent ±1σ for the ensemble over the 10-year period.

3.2.1 Optical depth

The ice extinction anomalies of G4-Base that are calculated in the ULAQ-CCM are negative in the whole UT (Fig. 12a, b) due to the decreasing number density of the particles caused by the reduced vertical velocities in the SG dynamical conditions (see Figs. 7–8). Although the UT cooling in G4 tends to partially offset the effects of the updraft decrease on the ice particle number density, the overall impact is of a general decrease of the UT ice extinction and is even more pronounced than in G4K where the tropospheric cooling is not taken into account. In the latter case, however, the particle effective radius is larger than in G4, as discussed above for Fig. 10. These size distribution changes affect not only ice extinction, but also the shortwave and longwave radiative responses per unit optical depth (see ahead Sect. 3.2.2).

Following the procedure described in Sect. 2.2 (see Eq. 3), an evaluation of the model calculated ice extinction profiles is attempted (Fig. 12c, d). This is made using indirectly derived values from the MERRA-2 and ERA5 ice mass mixing ratio and the ULAQ-CCM effective radius, as in Eq. (6) below. Here, χ_ext,i is the ice extinction at the ith vertical layer and ρ_atm,i is the atmospheric mass density at the same vertical layer:

$$\begin{array}{}\text{(6)}& {\mathit{\chi }}_{\mathrm{ext},i}=Q_{\mathrm{ext}}{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{atm},i}}{{\mathit{\rho }}_{\mathrm{ice}}}}{\displaystyle \frac{{\mathit{\chi }}_i}{r}}\end{array}$$

The ULAQ-CCM tropical underestimation of the ice extinction below 13 km is consistent with that of the ice number density and is partly justified by the specific assumptions made on cirrus cloud formation in the model, as pointed out in the discussion of Fig. 4.

The net result on the ice optical depth (i.e., the vertical integral of ice extinction) is shown in Fig. 13. In general, a latitude-dependent OD reduction comparable to that found in Kuebbeler et al. (2012) is present in G4K, while in the G4 case (as expected from the extinction anomalies) a further decrease is calculated mainly in the tropics, even though the UT temperatures are cooler. The effects regarding the temperature and updraft cannot be easily separated, but the colder tropospheric temperatures in G4 with respect to G4K reduce the particle size increase respect to the Base case, producing an additional decrease in the optical depth. The coupled effects of the velocity and temperature anomalies on the ice particle number density and size produce the most relevant impact in our study, pointing out the importance of allowing surface temperatures to respond to the stratospheric aerosol radiative forcing.

Figure 13 Zonally and time-averaged values of the ice optical depth (λ=0.55 µm) for the ULAQ-CCM Base, G4 and G4K experiments (solid black, red and blue lines, respectively) (a) and Base case comparison with the MERRA-2 and ERA5 indirectly derived values (dashed red and blue line) (b). The model results are for years 2030–2069; the MERRA-2 and ERA5 data are for years 2003–2012. The shaded area represents ±1σ for the ensemble over the 40-year period 2030–2069.

3.2.2 Consequences on radiative forcing

The well-tested radiative transfer code online in the ULAQ-CCM (Chou et al., 2001; Randles et al., 2013; SPARC, 2013) has been used to calculate the shortwave and longwave components of the tropopause radiative forcing due to SG aerosols (direct forcing) and to UT ice changes (indirect forcing). As discussed so far, the latter are largely produced by the SG-driven dynamical perturbations on the homogeneous freezing process for ice formation. The ice radiative effects have been calculated using up-to-date wavelength-dependent refractive index available in the literature (Warren, 1984; Warren and Brandt, 2008; Curtis et al., 2005) and compared against previous results under similar conditions, such as those by Schumann et al. (2012). All the radiative calculations shown in this section have been performed offline with the same radiative transfer code as the one present online in the ULAQ-CCM model, in order separate the effects of the single components analyzed.

The results are shown separately for the G4 and G4K experiments, both with respect to the RCP4.5 Base case. Following the previously discussed thinning of the UT ice clouds, a positive SW RF is calculated because of the decreased scattering of the incoming solar radiation by the ice particles. However, such an effect is largely covered by the negative LW RF due to a lessened capacity of the ice particles to trap outgoing planetary radiation; therefore, the obtained net effect on RF is negative, as shown in Table 3 and Fig. 14. This indirect negative RF is smaller but still significant when compared to the negative direct net RF due to the SG aerosols (∼30 % of it).

It is interesting to note that the shortwave component of the ice RF is indeed smaller than the longwave component, however, not as much as one could expect from the very different normalized RFs (i.e., forcing per unit OD) at a given particle radius. The reason is that both the SW and LW normalized RFs are decreasing with the increasing particle radius, but the relative changes of these normalized RF components are significantly different between the SW and LW. According to our radiative calculations, the SW normalized values decrease (in magnitude) from −12.1 to −5.7 W m⁻² (−53 %) with the ice effective radius increasing from 15 to 40 µm, whereas the instantaneous LW normalized RF remains quasi-constant at an average value of +53 W m⁻², with a smooth 3 % decrease over the same radius interval. The resulting SW RF is then controlled not only by the negative OD changes (−0.020 in G4 and −0.012 G4K) but also by the magnitude of the particle radius increase, which is larger in G4K than in G4, and larger in both perturbed cases compared with the Base case (see discussion of Fig. 10).

Table 3 Top three rows: globally and time-averaged values of the upper tropospheric ice optical depth changes and RF differences (W m⁻²) between the SG perturbed experiments and the RCP4.5 Base case due to changes in ice crystal concentration and size. Middle three rows: globally averaged values of stratospheric sulfate aerosol optical depth changes and RF differences (W m⁻²) defined as above but due to changes in aerosol concentration and size. Bottom three rows: total OD and RF changes (i.e., ice + sulfate). All results are for all-sky conditions (i.e., including the presence of background cloudiness) and with an 8 Tg-SO₂ yr⁻¹ injection. The RFs are calculated at the tropopause with temperature adjustment. The time is averaged over the years 2030–2069.

Figure 14 Zonally averaged cloud adjustments to the clear-sky SG aerosol RF (W m⁻²), as a function of latitude (time-averaged over the years 2030–2069). See legends for line meaning. The positive adjustment due to (passive) background clouds (black solid line for G4, black dashed line for G4K) shows the net value (SW+LW), which is largely controlled by the SW contribution (see Table 4 and Fig. S3).

Table 4 succinctly presents the globally and time-averaged ULAQ-CCM results for the cloud adjustments of the clear-sky RF components due to the SG stratospheric aerosols. The SW and LW cloud adjustments are roughly comparable to the those calculated in Kuebbeler et al. (2012) (+1.11 and −0.51 W m⁻², respectively, calculated at the top of atmosphere for an SG experiment with a 5 Tg-SO₂ yr⁻¹ injection). These numbers could be compared with those obtained in the ULAQ-CCM G4K case (although for an 8 Tg-SO₂ yr⁻¹ injection), i.e., +1.51 and −0.73 W m⁻² for SW and LW, respectively, with a net value of +0.52 W m⁻² compared with +0.60 W m⁻² in Kuebbeler et al. (2012).

In the (more realistic) G4 simulation performed by the ULAQ-CCM model, the SW cloud adjustment is only slightly smaller than in the G4K, while a significantly larger negative LW component is calculated. This results in a net adjustment of +0.52 W m⁻² in the G4 case compared with +0.72 W m⁻² in the G4K experiment. A latitude-dependent view of these results is presented in Fig. 16. The black solid line shows the net positive adjustment (SW+LW) due to the mere presence of background clouds, which substantially alter the radiative fluxes (see also Kuebbeler et al., 2012; Schulz et al., 2006; Stier et al., 2013). These clouds are kept fixed in the ULAQ-CCM model, using climatological values, and thus do not present changes under the G4 scenario. An estimate of the all-sky RF contribution due to SG-driven changes of background clouds is beyond the scope of this study. According to our model calculations, the negative LW is the dominant component of the cloud adjustment due to cirrus ice thinning, and this is particularly true for the more realistic G4 simulation. In this latter case, significantly larger values of the LW adjustment are found over the tropics with respect to G4K, which is consistent with the ice extinction profile changes in Fig. 12a.

Table 4 Rearrangement of the results presented in Table 3, with the calculated cloud adjustments (bottom three rows) to clear-sky RF components (top three rows). The cloud adjustments for the SW and LW RF contributions are shown separately for the presence of background atmospheric clouds (left) and for cirrus thinning (right): the former is calculated as the difference between the all-sky and clear-sky aerosol RFs, with the all-sky including the background warm clouds and fixed UT ice clouds.

Further information regarding the model calculated RFs is presented in Fig. S3, where we show both the clear-sky latitudinal distribution of the sulfate aerosol RF (Fig. S3a) for G4 and G4K and the LW and SW cloud adjustment due to the presence of background clouds for G4 and G4K (Fig. S3b).

4 Conclusions

Sulfate geoengineering is currently considered to be one of the most promising solar radiation management (SRM) techniques. One reason for this (and unlike other methods) is that past explosive volcanic eruptions in the tropical belt are a natural proxy for stratospheric sulfate injection. However, this does not mean that SG does not still pose some scientific questions that need to be thoroughly answered, as has been pointed out by MacMartin et al. (2016). For instance, models still show many significant differences regarding the confinement of stratospheric sulfate aerosols in the tropical pipe (Pitari et al., 2014).

In recent years, some experiments have been proposed in which SG is used to meet different climate targets (MacMartin et al., 2017; Kravitz et al., 2017). However, to properly utilize SG, a clear understanding is needed of how multiple side effects of this technique can modify the net RF (Visioni et al., 2017a). While some of these effects produce a negligible difference in forcing, such as those from gas species perturbations (CH₄, O₃, stratospheric H₂O) (Visioni et al., 2017b), this might not be the case for changes produced in the formation of thin cirrus ice clouds.

This latter indirect effect has already been analyzed in two previous works. Cirisan et al. (2013) looked at the potential impact of IN changes in the UT, finding a negligible positive TOA (top of the atmosphere) forcing (+0.02, up to 0.04 W m⁻²) due to the number density increase of H₂SO₄–H₂O aerosols transported down to the UT from the lower stratosphere. Kuebbeler et al. (2012), in comparison, studied the effects of dynamical changes caused by the aerosol-induced stratospheric warming and their consequences on UT ice formation via homogeneous freezing. They found a considerable negative TOA forcing in the longwave spectrum (−0.51 W m⁻²), which was greatly attributable to the SG-induced ice optical depth reduction. In the present study, we focus on these same indirect dynamical effects, adding the potential impact of the SG aerosol-induced surface cooling (G4 experiment), which was not explicitly considered in Kuebbeler et al. (2012). Their approach was also included for comparison in our study, by means of a sensitivity study (G4K) conducted using the ULAQ-CCM model, where surface temperature remained fixed at the RCP4.5 baseline values so that we could more precisely quantify the surface cooling impact on the UT thin cirrus clouds.

Figure 15 Schematic summary of the sulfate geoengineering impact on the dynamical processes driving changes of upper tropospheric ice particle formation through homogeneous freezing.

Figure 16 Illustration of the sulfate geoengineering impact on cirrus ice particles formed through freezing and schematic representation of ice and aerosol changes in radiative fluxes.

A compact view of the SG effects on UT ice formation is presented in Fig. 15. On one hand, the aerosol-induced stratospheric warming and surface cooling combined produce a further atmospheric stabilization with an even larger reduction in tropospheric updraft compared with the G4K case. This lowers the UT probability of ice supersaturation, with less favorable conditions especially for homogeneous freezing. On the other hand, this ice formation limiting effect is partially counterbalanced by the convectively driven tropospheric cooling, which is not observed in the G4K case.

The resulting changes in ice particle number density and size distribution, when combined, translate into a globally averaged decrease of the ice optical depth ($\mathrm{\Delta }\mathit{\tau }=-\mathrm{0.024}$, at λ=0.55 µm), i.e., −6 % of the baseline OD. This reduction is larger than the one in G4K relative to the Base case ($\mathrm{\Delta }\mathit{\tau }=-\mathrm{0.012}$, −3 %), which points to the dominant and controlling role of the reduced updraft velocities. According to our model results, these OD changes (coupled to increases in ice particle effective radii) translate to net tropopause RFs of −0.29 and −0.14 W m⁻², for G4 and G4K experiments, respectively, produced only by the cirrus ice thinning effect of SG. These two cloud adjustments result from a combination of the SW and LW RF contributions, which account for +0.50 and +0.35 W m⁻² in the SW (for G4 and G4K, respectively) and −0.79 and −0.49 W m⁻² in the LW (again for the respective G4 and G4K experiments).

We can compare these ice thinning forcing contributions with the net tropopause all-sky RF produced by the stratospheric SG aerosols, i.e., of −1.17 and −1.24 W m⁻², for the G4 and G4K experiments, respectively. According to our model, the net negative RF due to the cirrus ice cloud thinning is close to 25 % (in G4) of the direct effect of the sulfate particles themselves. This might have consequences regarding the definition of the sulfate injection efficiency in terms of RF per Tg-S yr⁻¹ injected, especially if such efficiency is used to determine the amount of SO₂ that needs to be injected into the stratosphere to achieve climate targets (MacMartin et al., 2017; Kravitz et al., 2017).

Figure 16 summarizes the thermodynamical processes leading to the changes in cirrus ice formation and the radiative response caused by these changes in the Earth's radiative balance, as analyzed in detail in this paper, in addition to the direct radiative effect of the sulfate particles.

Furthermore, one last consideration is necessary regarding the RFs in the SG scenarios and the unperturbed atmosphere, more specifically, with respect to the cloud adjustment to clear-sky RFs due to the stratospheric sulfate aerosols. In our fully interactive aerosol simulation (G4), we obtain a total cloud adjustment (from both cirrus ice thinning and passive background clouds) of +0.60 W m⁻² due to compensating for large adjustments in the LW and SW. The SW adjustment results in part from the mere presence of (passive) background clouds and in part from the changing size distribution of UT ice particles. The increasing particle size is more pronounced in the partially interactive aerosol simulation (G4K), and thus produces a larger positive SW contribution (+0.78 W m⁻²).

This latter value is comparable to that calculated in the similar experiment of Kuebbeler et al. (2012) (+0.60 W m⁻², with a 5 Tg-SO₂ injection). This means that the lower stratospheric warming produced by the SG aerosols acts indirectly on atmospheric dynamics with a strong feedback on the UT cirrus clouds; therefore, a simple reduction of the incoming solar radiation is not a good proxy for the eventual injection of sulfate particles into the stratosphere. When the aerosol-induced surface cooling is coupled to the lower stratospheric warming, the net cloud adjustment is significantly reduced; however, the clear-sky balance of the SW and LW RF contributions is greatly altered by the presence of background clouds coupled to the UT ice thinning.

One important caveat to the conclusions of this study, is that the physical processes behind the UT ice particle formation are highly idealized in our parameterization. Nonetheless, the results it produces in the reference (historical) simulation are generally comparable with the MERRA-2 and ERA5 reanalysis and some satellite data. In addition, the calculated SG dynamical anomalies in the stratosphere are consistent with those from other modeling studies (Pitari et al., 2014; Niemeier and Schmidt, 2017). Finally, considering the consistency with the findings from the study of Kuebbeler et al. (2012), we may reasonably conclude that our results regarding the thinning of the UT ice clouds under SG conditions are sufficiently robust. However, considering how complex the balance between the UT ice formation changes and their radiative forcing is (Mitchell et al., 2008), the results from this study cannot be considered conclusive and exhaustive. Additional results using different and more complete physical parameterizations (both regarding the ice formation processes and a wider range of updraft velocities), and an online ocean coupling, may help clarify the net contribution of ice clouds in a sulfate geoengineering scenario.