The effect of ash, water vapor, and heterogeneous chemistry on the evolution of a Pinatubo-size volcanic cloud

We employ the ECHAM5/MESSy2 atmospheric chemistry general circulation model (EMAC) that incorporates calculations of gas-phase and heterogeneous chemistry coupled with the ozone cycle and aerosol formation, transport, and microphysics to calculate the 1991 Pinatubo volcanic cloud. We considered simultaneous injections of SO₂, volcanic ash, and water vapor. We conducted multiple ensemble simulations with different injection configurations to test the evolution of SO₂, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, ash masses, stratospheric aerosol optical depth, surface area density (SAD), and the stratospheric temperature response against available observations. We found that the volcanic cloud evolution is sensitive to the altitude where volcanic debris is initially injected and the initial concentrations of the eruption products that affect radiative heating and lofting of the volcanic cloud. The numerical experiments with the injection of 12 Mt SO₂, 75 Mt of volcanic ash, and 150 Mt of water vapor at 20 km show the best agreement with the observation aerosol optical depth and stratospheric temperature response. Volcanic water injected by eruptive jet and/or intruding through the tropopause accelerates SO₂ oxidation. But the mass of volcanic water retained in the stratosphere is controlled by the stratospheric temperature at the injection level. For example, if volcanic materials are released in the cold point above the tropical tropopause, most of the injected water freezes and sediments as ice crystals. The water vapor directly injected into the volcanic cloud increases the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass and stratospheric aerosol optical depth by about 5 %. The coarse ash comprises 98 % of the ash injected mass. It sediments within a few days, but aged submicron ash could stay in the stratosphere for a few months providing SAD for heterogeneous chemistry. The presence of ash accelerates the SO₂ oxidation by 10 %–20 % due to heterogeneous chemistry, radiative heating, lofting, and faster dispersion of volcanic debris. Ash aging affects its lifetime and optical properties, almost doubling the ash radiative heating. The 2.5-year simulations show that the stratospheric temperature anomalies forced by radiative heating of volcanic debris in our experiments with the 20 km injection height agree well with observations and reanalysis data. This indicates that the model captures the long-term evolution and climate effect of the Pinatubo volcanic cloud. The volcanic cloud's initial lofting, facilitated by ash particles' radiative heating, controls the oxidation rate of SO₂. Ash accelerates the formation of the sulfate layer in the first 2 months after the eruption. We also found that the interactive calculations of OH and heterogeneous chemistry increase the volcanic cloud sensitivity to water vapor and ash injections. All those factors must be accounted for in modeling the impact of large-scale volcanic injections on climate and stratospheric chemistry.

1 Introduction

Strong explosive volcanic eruptions are the major natural cause of climate variability on both global and regional scales (Robock, 2000). Volcanic explosions inject a mixture of SO₂, volcanic ash, water vapor, halogens, and other tracers into the stratosphere. The injected volcanic materials scatter and absorb incoming solar and outgoing terrestrial radiation, warming the stratosphere and cooling the Earth's surface and the lower troposphere (Hansen et al., 1992; Stenchikov et al., 1998; Kirchner et al., 1999; Robock, 2000; Soden, 2002; Shindell et al., 2001). Stratospheric warming and tropospheric cooling caused by the radiative impact of volcanic aerosols yield to changes in atmospheric circulation (Stenchikov et al., 1998; Kirchner et al., 1999; Ramachandran et al., 2000), affect El Niño–Southern Oscillation (ENSO) (Predybaylo et al., 2017), and force a positive phase of the Arctic Oscillation (AO) (Graf et al., 1993; Kodera and Kuroda, 2000; Mao and Robock, 1998; Kodera and Kuroda, 2000; Stenchikov, 2002; Shindell, 2004; Stenchikov et al., 2006; Karpechko et al., 2010) causing boreal winter warming in middle and high latitudes over Eurasia and North America (Stenchikov et al., 2004; Thomas et al., 2009; Poberaj et al., 2011).

Recent research shows that the mechanisms of stratosphere–troposphere dynamic interaction are complex, and large simulation ensembles are required to detect and attribute dynamic responses reliably. Stenchikov et al. (2006), Driscoll et al. (2012), and Charlton-Perez et al. (2013) showed that the Intergovernmental Panel on Climate Change (IPCC) models have the problem of producing a stronger northern polar vortex in response to low-latitude volcanic eruptions. Conveying this signal to the surface is even more problematic. For example, Polvani et al. (2019) concluded that the positive AO phase after the Pinatubo eruption appeared only by chance. Toohey et al. (2014) further elaborated on the planetary-wave-based mechanism of winter warming after large low-latitude eruptions. Bittner et al. (2016) and Azoulay et al. (2021) showed that a stronger eruption could more reliably force a positive phase of the AO.

The eruption of Mount Pinatubo in the Philippines on 15 June 1991 had an explosivity index of VEI = 6 and caused the largest climate impact in the 20th century. It was also the largest eruption that affected a densely populated area. The observed global mean visible optical depth from the Pinatubo eruption reached 0.15. It was about 2 times higher than the second largest eruption in the 20th century: El Chichón in 1982 (Dutton and Christy, 1992). The 1991 Mt. Pinatubo eruption is also the best-observed explosive event with a detected significant climate impact. It has been documented by satellite instruments (McCormick, 1987; Long and Stowe, 1994), ground-based lidars and sun photometers (Antuna et al., 2002, 2003; Good and Pyle, 2004; Nagai et al., 2010; Dutton and Christy, 1992; Thomason, 1992) , and airborne aerosol counters (McCormick et al., 1995; Pueschel et al., 1994; Borrmann et al., 1995; Deshler, 2003). Mount Pinatubo produced about 5 km³ of magma. According to observations, three main volcanic explosions on 15 June 1991 spread volcanic ash and gases over an area of 300 000 km²

The SO₂ mass emitted by the Mount Pinatubo eruption was estimated using the Stratospheric Aerosol and Gas Experiment (SAGE), the TIROS Operational Vertical Sounder (TOVS), and the Total Ozone Mapping Spectrometer (TOMS) retrievals (Guo et al., 2004a; Rose et al., 2006; Sheng et al., 2015; Krueger et al., 1995; Fisher et al., 2019). In addition to SO₂, Pinatubo injected tens of megatons of water vapor and volcanic ash into the stratosphere (Guo et al., 2004a; Nedoluha et al., 1998; Joshi and Jones, 2009).

In the stratosphere, SO₂ is oxidized by the OH radical to form sulfuric acid, which then binary-nucleates in the presence of water to form sulfate aerosol. The primary source of OH in the stratosphere is ozone photolysis by ultraviolet radiation. This reaction forms oxygen and atomic oxygen in the excitation state (O^1d), which interacts with water vapor to form OH radicals. Thus, the SO₂ oxidation is controlled by the abundance of OH, which depends on the concentration of stratospheric water vapor (Lovejoy et al., 1996). The co-injection of water vapor with SO₂ accelerates the formation of sulfuric acid (LeGrande et al., 2016). The online calculation of OH is essential to correctly reproduce the dynamics of sulfate aerosol mass (Clyne et al., 2021; Stenchikov, 2021), and this has been neglected in some previous studies (Marshall et al., 2018; Niemeier et al., 2009; Oman et al., 2006).

The sulfuric acid resulting from SO₂ oxidation nucleates to form long-lived submicron sulfate droplets interacting with solar and terrestrial radiation. The radiative effect and lifetime of sulfate aerosols depend on their size distribution, which is not definitively established. Therefore, different Pinatubo modeling studies report a wide range of visible (0.5–0.6 µm) stratospheric aerosol optical depth (SAOD) for the same amount of injected SO₂. Brühl et al. (2015) obtained equatorial average SAOD = 0.38 compared to SAOD = 0.11 reported by LeGrande et al. (2016) for 17 Mt of injected SO₂. Niemeier et al. (2021) and Stenchikov et al. (2021) obtained similar SAOD which is consistent with observations for 17 Mt of injected SO₂. Dhomse et al. (2014), using a detailed aerosol microphysics model, found that in simulations of a Pinatubo-like eruption with 10 Mt of SO₂ injection, SAOD matches observations better than that with larger SO₂ emission. Mills et al. (2016) also reported that in their model a 10 Mt SO₂ injection produces the best fit to Pinatubo observations, while Sheng et al. (2015) and Sukhodolov et al. (2018) found that SAOD in their experiments with the emission of the 14 Mt of SO₂ best fits SAGE observations. Timmreck et al. (2018) conducted ensembles of simulations with perturbed parameters, including the mass of injected SO₂ and the injection height of volcanic debris, to quantify the uncertainties in the radiative forcing of the 1991 Mt. Pinatubo eruption.

Volcanic ash (tephra) comprises silicate and volcanic glass with traces of gas bubbles (Kremser et al., 2016). Ash particles have a wide range of sizes from sub-micrometers to millimeters (Rose and Durant, 2009) and highly irregular shapes. Large ash particles with radii r > 1 µm sediment relatively quickly (Niemeier et al., 2021, 2009; Stenchikov et al., 2021) and are believed to contribute little in the long-term evolution of a volcanic cloud. Fine ash particles with r < 1 µm disperse over vast distances and can survive in the stratosphere for several months (Pueschel et al., 1994; Zhu et al., 2020; Russell et al., 1996; Vernier et al., 2016). Still, their radiative effect is small because of their relatively small mass. However, Stenchikov (2021) showed that although most ash mass sediments during the first week after an eruption, ash solar and infrared (IR) heating and chemical/microphysical interactions with sulfate particles could affect the volcanic cloud formation and its long-term evolution. Ash particles could be coated by sulfate, becoming chemically aged (Muser et al., 2020; Zhu et al., 2020). They also uptake SO₂, thereby decreasing its abundance (Zhu et al., 2020). The coating and aging of ash particles increase their size, alter their optical properties, and increase their deposition velocities (Muser et al., 2020; Zhu et al., 2020). The enhanced ash sedimentation removes a portion of sulfate depending on the aging level. At the same time, stratospheric aerosol particles (ash and sulfate) provide surfaces for heterogeneous chemical reactions affecting stratospheric chemical composition (Muthers et al., 2015). Aerosol particles from volcanic eruptions increase the surface area density (SAD) and hence the rate of heterogeneous reactions involving ClONO₂ and N₂O₅. This damps the NO_x mixing ratios altering the OH stratospheric budget, which affects the rate of SO₂ oxidation (Prather, 1992; Kilian et al., 2020).

Volcanic ash absorbs and scatters solar shortwave (SW) and terrestrial longwave (LW) radiation. This significantly impacts the dynamics and radiation budget of the atmosphere in the first few days after an eruption, causing rapid lifting of volcanic debris (Stenchikov et al., 2021; Niemeier et al., 2009). The models that calculate the evolution of volcanic clouds and their impact on climate assume various SO₂ injection heights and initial plume composition (ash and water are often not injected), use different spatial–temporal resolutions, and treat ash–sulfate chemical interaction differently. The differences in physics and chemistry translate into the differences in volcanic cloud evolution and radiative effect.

Along with the SO₂ mass, the injection height of volcanic debris is a critical parameter for correctly simulating the dispersion of a volcanic cloud, as it is associated with the wind field that transports the volcanic plume. Different modeling studies assume different injection heights. The results show that the oxidation rate of SO₂ depends strongly on the injection height according to the availability of water vapor and OH radicals. Sheng et al. (2015) performed a sensitivity study for the initial mass and altitude of the injected SO₂ for the Pinatubo eruption and showed that a mass of 17 Mt of SO₂ gives the best agreement with the SAGE optical depth within a peak of the volcanic cloud between 18–21 km. The transient equilibrium height of the volcanic plume depends not only on the height of the initial injection but also on internal feedback mechanisms. Stenchikov et al. (2021) demonstrated that radiative heating by ash was lifting volcanic debris by 1 km d⁻¹ during the first week following the 1991 Pinatubo eruption. Muser et al. (2020) reported the lifting of a volcanic plume of the Raikoke eruption, 1 order of magnitude smaller than the 1991 Pinatubo eruption. Volcanic debris injections cause significantly weaker lofting when their initial distribution is quasi-zonal, as in Brühl et al. (2015), than when it is localized, as in LeGrande et al. (2016). In addition to processes considered in Stenchikov et al. (2021), we here explicitly calculate ash chemical aging, stratospheric ozone chemistry, and aerosol microphysical processes, accounting for the hygroscopic growth of sulfate/ash particles. However, we do not account for the heating by SO₂ because for the Pinatubo eruption, it is much weaker than radiative heating from ash and sulfate aerosols (Osipov et al., 2020, 2021; Stenchikov et al., 2021).

The underlying dynamic and chemical mechanisms of the large sensitivity of SAOD to the injection height have not been tested yet in a fully interactive model. The effects of injected “volcanic” water and chemical aging of volcanic ash on SO₂ oxidation rate and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ removal are not studied within the models with comprehensive gaseous and heterogeneous chemistry and detailed microphysics. Here we use the EMAC model with multi-phase chemistry along with detailed aerosol microphysics to study the evolution of a Pinatubo-size volcanic cloud. We account for the entire range of dynamic, chemical, and microphysical complexity of the processes governing the development of volcanic cloud to address the following science questions:

• How do the initial spatial distribution and height of injected volcanic debris affect the evolution of a volcanic cloud?

• What is the effect of heterogeneous chemistry on the SO₂ oxidation rate within a volcanic cloud?

• How does the co-injection of SO₂, water vapor, and ash affect volcanic cloud evolution?

• How does the aging of co-injected ash affect volcanic cloud development?

2 Data

To constrain the simulations and evaluate the model results, we use the SAGE data set with partially filled gaps compiled by Stratospheric Processes And their Role in Climate (SPARC) and published in the Assessment of Stratospheric Aerosol Properties (ASAP) report (Thomason and Peter, 2006). This data set provides the aerosol effective radius and aerosol extinction in UV, visible (0.525 µm), and near-IR (1.02 µm) wavelengths. The SAGE/ASAP SAOD is a zonal mean and assembled on a monthly basis. It is available from 70^∘ S to 70^∘ N with a 5^∘ resolution in latitude from 1984 to 1999. We refer to these data as SAGE/ASAP SAOD or R_eff. The SAGE observations of aerosol extinction contain multiple gaps in the initial stage of the volcanic cloud evolution because of the instrument's saturation (Thomason, 1992). The near-IR observations are of better quality than in the visible or ultraviolet (UV) range (Stenchikov et al., 1998). Therefore, to obtain the visible SAGE/ASAP SAOD, we scale near-IR SAOD using the Ångström exponent obtained from our simulations, similar to Stenchikov et al. (1998). We also use the Advanced Very High-Resolution Radiometer (AVHRR) SAOD at 0.63 µm (Long and Stowe, 1994). The AVHRR observations are collected over the oceans at 0.1^∘ × 0.1^∘ horizontal resolution for cloud-free conditions in daylight time. The AVHRR aerosol optical depth (AOD) is measured for the entire atmospheric column, including the troposphere. To obtain the AVHRR stratospheric AOD, we calculate the AVHRR AOD monthly climatology for the pre-Pinatubo period of 1985–1990 and subtract it from the total AOD for the Pinatubo period. Unfortunately, this can introduce uncertainty due to the high variability of tropospheric AOD. We refer to the visible SAOD obtained from the SAGE near-IR observations as the scaled SAGE/ASAP SAOD.

Krueger et al. (1995) estimated the mass of SO₂ during the first 15 d after the 1991 Pinatubo eruption based on the Total Ozone Mapping Spectrometer (TOMS) observations. They concluded that the mass of the initially emitted SO₂ was 15 ± 3 Mt. Guo et al. (2004a) later estimated the emitted mass of SO₂ to be 14–20 Mt. Recent estimates reduce the initial SO₂ mass to 12 Mt (Fisher et al., 2019). Estimates of SO₂ mass using retrievals from the TIROS Optical Vertical Sounder/High-Resolution Infrared Radiation Sounder/2 (TOVS) on the Television Infrared Observation Satellite (TIROS) suggest that the initial SO₂ mass was 19 ± 4 Mt (Guo et al., 2004a). However, the TOVS retrievals are less accurate than TOMS because they are affected by sulfate aerosol absorption in the IR range. ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass has also been estimated using the High-Resolution Infrared Radiation Sounder/2 (HIRS/2) (Guo et al., 2004a). The estimated sulfate aerosol mass depends on the aerosol size distribution, which is not well known, and this adds uncertainties to the estimated ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass.

Volcanic ash mass was estimated for the first few days after the Pinatubo eruption using the AVHRR and HIRS/2 observations (Guo et al., 2004b). HIRS/2 detected 80 Mt of fine ash in the atmosphere on the first day after the eruption. AVHRR ash retrievals evaluate the spectral contrast of radiance (aerosol index) to distinguish between absorbing aerosols, such as volcanic ash, and non-absorbing aerosols, such as sulfate. However, the retrieval algorithm does not consider particles smaller than 1 µm (Guo et al., 2004b).

We obtain the stratospheric temperature response to the 1991 Pinatubo eruption from the Modern-Era Retrospective analysis for Research and Applications (MERRA2) reanalysis data available on a 0.5^∘ × 0.625^∘ horizontal grid and 72 vertical levels from the surface to 0.01 hPa (Gelaro et al., 2017). To reproduce the effect of the 1991 Pinatubo eruption, MERRA2 assimilates observations from different satellite sensors such as TOVS and the Meteosat Second Generation (MSG) satellite as well as the Microwave Limb Sounder (MLS). The MERRA2 temperature fields are consistent with the observations reported by Labitzke and McCormick (1992). The reanalysis intercomparison S-RIR (Fujiwara et al., 2022) shows that the MERRA2 stratospheric temperature anomalies caused by the 1991 Pinatubo injection resemble the observations well, although the absolute stratospheric temperature is slightly underestimated. In this study, the MERRA2 temperature anomalies were calculated with respect to 1985–1990 climatology.

3 Model

Here we employ the ECHAM5/MESSy2 atmospheric chemistry model, EMAC (Joeckel et al., 2005, 2006, 2010). EMAC is a modular model based on submodels that describe processes in the stratosphere, the middle atmosphere, and the troposphere, accounting for anthropogenic emissions and interactions with oceans and land (Joeckel et al., 2010). EMAC has been used to study impacts of volcanic stratospheric aerosols on climate and stratospheric circulation (Brühl et al., 2012, 2015; Bingen et al., 2017; Löffler et al., 2016; Kilian et al., 2020) as well as dust aging and dust–air pollution interactions in the troposphere (Abdelkader et al., 2015, 2017; Klingmüller et al., 2019, 2020).

Klingmüller et al. (2014)Pozzer et al. (2006)Roeckner et al. (2006)Tost et al. (2006b, 2010)Tost et al. (2006b)Kerkweg et al. (2006a)Pringle et al. (2010)Landgraf and Crutzen (1998)Sander et al. (2014, 2019)Tost et al. (2007)Sander et al. (2011)Kerkweg et al. (2006b)Kerkweg et al. (2006b)Roeckner et al. (2006)Joeckel et al. (2006)Tost et al. (2006a)Kerkweg et al. (2006a)Kerkweg et al. (2006b)Joeckel et al. (2006)

Table 1 List of EMAC submodels used in this study. A complete list of all EMAC submodels can be found in Joeckel et al. (2010).

The Modular Earth Submodel System (MESSy) links the various submodels. The submodels comprise AEROPT, CLOUD, CONVECT, CVTRANS, DDEP, GMXE, JVAL, LNOX, MECCA, OFFEMIS, ONEMIS, RAD4ALL, SCAV, SEDI, TNUDGE, and TROPOP. Table 1 shows the submodels used in this study indicating their functionality, while the detailed description of all EMAC submodels can be found in Joeckel et al. (2010). We configure EMAC using MESSy version 2.52 with the fifth-generation European Centre Hamburg Atmospheric general circulation Model version 5.3.02, ECHAM5 (Roeckner et al., 2006), and employ the same chemistry and aerosol microphysics setup as in Brühl et al. (2012, 2015). For vertical approximation, we use 90 sigma hybrid levels from the Earth surface up to 0.01 hPa and a T42 spectral approximation horizontally corresponding to 2.8^∘ grid spacing at the Equator both in longitude and latitude. Varying monthly sea surface temperature and sea ice are prescribed from the AMIPII data set (Taylor et al., 2000).

We apply the quasi-biennial oscillation (QBO) submodel to capture the observed phase of QBO and account for its effect on the stratospheric circulation, similar to Stenchikov et al. (2004). No other constraints are imposed on the model dynamics.

The emission inventory comprises the sources of greenhouse gases, NO_x, CO, non-methane volatile organic compounds (NMVOSs), NH₃, SO₂, black carbon (BC), and organic carbon (OC). We also account for the background emission of CFCs, halogens, and halons. The emissions are monthly means and geographically distributed according to the 2009 Emissions Database for Global Atmospheric Research version 4 (EDGAR4) emission inventory and the Global Fire Emissions Database (GFED) version 3 (van der Werf et al., 2010). We account for the dimethylsulfide (DMS) and carbonyl sulfide (OCS) emissions similarly to Brühl et al. (2015). We employ 230 gas-phase chemical reactions, 76 photolytic reactions, and 12 heterogeneous reactions for 159 species to calculate atmospheric composition. The photolysis rates are calculated within the model for the spectral range of 178.6 nm ≤ λ ≤ 752.5 nm accounting for gaseous absorption (O₃ and O₂), Rayleigh scattering, absorption, and scattering by aerosols and clouds (Landgraf and Crutzen, 1998; Sander et al., 2011). In this setup, the photolysis rates are not coupled to volcanic aerosol. The model calculates the instantaneous radiative forcing using double radiation calls, with and without aerosols. Aerosol microphysics and chemistry are called every model time step, while the radiation submodel is called every third time step.

3.1 Stratospheric sulfate chemistry

Volcanic sulfate results from the oxidation of SO₂ by OH in the presence of water vapor. OH is produced by ozone photolysis by UV radiation with wavelengths of less than 0.242 µm. This reaction forms O₂ and excited oxygen O(^1d). The excited oxygen radical interacts with water to form the hydroxyl radical OH, which oxidizes SO₂ in two steps to form sulfate. At the first step, OH oxidizes SO₂ to form SO₃ and HO₂. At the second step, SO₃ interacts with water molecules to form sulfuric acid (Seinfeld and Pandis, 2006). The rate of reaction depends on the concentration of water molecules that are also the reactants (Burkholder et al., 2015). Therefore, higher water vapor concentration significantly increases the formation rate of sulfuric acid. New sulfate particles are generated by the binary nucleation of sulfuric acid and water molecules. Thus, the formation of sulfate particles in a volcanic cloud depends strongly on water vapor concentration. The models that do not parameterize nucleation explicitly are less sensitive to the abundance of water vapor in a volcanic cloud than those that do (LeGrande et al., 2016).

Figure 1 Schematic representation of the initial size distribution of aerosol modes in EMAC (nucleation, Aitken, accumulation, and coarse). The threshold radii R1, R2, R3, and R4 control particle exchange between the modes. Initially 1.5 Mt of volcanic ash was injected in accumulation mode and 73.5 Mt in coarse mode.

3.2 Aerosol microphysics

The aerosol setup in EMAC has been described in detail in Pringle et al. (2010), Tost et al. (2010), de Meij et al. (2012), Pozzer et al. (2012), Brühl et al. (2015), and Abdelkader et al. (2015, 2017). We use the aerosol microphysics submodel GMXe (Pringle et al., 2010), coupled to the gas–aerosol partitioning scheme ISORROPIA-II (Fountoukis and Nenes, 2007) and heterogeneous chemistry schemes (Sander et al., 2005, 2011, 2019). Aerosol size distributions in the model are approximated by seven lognormal modes: four soluble modes (nucleation, Aitken, accumulation, coarse) and three insoluble modes (Aitken, accumulation, coarse). In our simulations, sulfate particles are represented by soluble modes, while ash is initially considered insoluble until it ages; i.e., five monolayers of sulfate particles coat the ash particle. The modes' median radii change in time during aerosol microphysical transformations, but the widths of the modes remain fixed. The median radii for the three insoluble modes and dry cores of the four soluble modes initially are equal to 0.0015, 0.025, 0.25, and 2.5 µm for nucleation, Aitken, accumulation, and coarse modes, respectively. The geometric standard deviations of the lognormal distributions for the above modes are 1.59, 1.59, 1.49, and 1.70, respectively (Brühl et al., 2015).

Aerosols in soluble modes evolve by uptake or loss of water and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ molecules and coagulation. The hygroscopic growth of ash is only allowed in a soluble mode (Abdelkader et al., 2015). The mass of large or fine aerosol particles in the size distribution tails is assigned to a corresponding neighboring mode when the mode's median radius reaches a certain threshold. The aerosol modes and the thresholds are schematically shown in Fig. 1. In our simulations, we choose threshold radii equal to 0.0005, 0.006, 0.07, and 1.6 µm for the nucleation, Aitken, accumulation, and coarse modes, respectively, as in Brühl et al. (2015).

3.3 Volcanic ash – model implementation

We introduced a new “ash” tracer to account for volcanic ash in the EMAC setting. We assume the ash density to be 2400 kg m⁻³. Ash is similar to mineral dust, comprising mainly silicates (SiO₂). Therefore, for calculating chemical aging, we assume that ash particles have the same water uptake and accommodation coefficients as dust particles (Abdelkader et al., 2015).

In our simulations, high-density ash particles sediment faster than pumice assumed in Zhu et al. (2020). Zhu et al. (2020) considered that the Kelud eruption emitted 100 times less volcanic material than the 1991 Pinatubo eruption; therefore, even the long-lived pumice ash produced negligible radiative heating. Stenchikov et al. (2021) showed that applying the assumption about long-lived pumice ash for the larger volcanic explosions like the 1991 Pinatubo eruption could cause unrealistic overheating of the stratosphere. Volcanic ash is removed from the stratosphere mainly by gravitational sedimentation. Sedimentation parameterization in EMAC utilizes the Walcek scheme (Walcek, 2000; Kerkweg et al., 2006a). Ash scavenging in the troposphere is implemented in EMAC by Tost et al. (2006a) and is fully coupled with aerosol and gas-phase chemistry.

For a complete representation of chemical aging, we use a comprehensive chemistry scheme that enables the production of the primary inorganic acids which contribute to the chemical aging of ash particles (Metzger et al., 2016).

Ash particles scatter and absorb solar and terrestrial radiation. To calculate their optical properties we choose the volcanic ash's complex refractive index according to Pollack et al. (1973), Vogel et al. (2017), and Stenchikov et al. (2021). For visible wavelengths, the ash refractive index RI = 1.53 + 0.004i. Ash is more absorbing at UV, near-IR, and IR wavelengths than at visible ones. Table S1 in the Supplement shows the volcanic ash refractive index as a function of wavelength. Figure S8 shows the refractive indices used in the EMAC model for different aerosols as a function of wavelength.

3.4 Aerosol radiative effect

We use the AEROPT submodel to calculate extinction, single-scattering albedo, and asymmetry parameter, the aerosol optical properties required for the radiative transfer calculations. It is assumed that different types of aerosols are mixed internally so that the refractive index of the mixture is calculated from the volume fractions of the aerosol components. The sensitivity to this assumption is discussed in detail by Klingmüller et al. (2014). The optical properties are calculated for each aerosol mode independently. The RAD submodel calculates radiative transfer (Roeckner and Coauthors, 2003). The Fouquart and Bonnel (1980) scheme is used for calculating shortwave radiation, while longwave radiation is calculated using a rapid radiative transfer model (RRTM) (Iacono et al., 2008). Scattering of the IR radiation by aerosols is neglected. RAD accounts for shortwave and longwave absorption of water vapor, clouds, O₃, CH₄, N₂O, CO₂, CFCs, and aerosols. Table 2 shows the shortwave and longwave bands used in the radiative transfer calculations in EMAC. For comparison with observations, we consider the first two SW bands in Table 2 as visible and near-IR.

Table 2 Shortwave and longwave bands used in the radiation transfer calculations.

Table 3 Description of experiments. The experiments are labeled according to the initial injection size and constituents of the injected plume. All experiment, the “1 × 1” format represents injection in one grid box, 3s10 represents zonal injection with 10 grid points in the latitude direction, the letter “s” denotes that only SO₂ was injected (dry injection), the letter “w” denotes that SO₂ and water vapor were injected (wet injection), and the va0 injection of volcanic ash with no aging and va1 is the aging case. For the 1w1,va0, and va1 experiments, 15 and 150 Mt of water vapor are injected. Each has five ensemble members.

¹ One grid box: 280 km × 280 km with a thickness of 1 km at 17 and 20 km altitude and 0.5 km at 25 km altitude. ² Ten grid boxes in latitude and 10 grid boxes in height (5km thickness).

Aquila et al. (2012)English et al. (2013)Bânda et al. (2013)Dhomse et al. (2014)Bândă et al. (2015)Sheng et al. (2015)Mills et al. (2016)

Table 4 List of the studies that simulated interactive chemistry for the Pinatubo case and the injected SO₂ height, maximum AOD, and the time (in months) for the maximum AOD.

NA: not available.

4 Experimental setup

Table 3 shows the complete set of numerical experiments conducted in this study. The control experiment (ctrl) describes the state of the atmosphere from 1990 to 2000, assuming there were no volcanic eruptions. In the perturbed simulations (with the injection of volcanic aerosols), we emit 17 Mt of sulfur dioxide as in Stenchikov et al. (2021), except in one experiment, specifically marked and used to study sensitivity to SO₂ emission mass. The perturbed simulations were calculated from 1 June 1991 to 31 December 1994 with a 1-year spin-up not included in the analysis.

Along with SO₂, we consider co-injections of water vapor and ash. For ash, we adopt the same initial size distribution as in Niemeier et al. (2009) and Stenchikov et al. (2021). We redistribute the total emitted fine ash mass of 75 Mt (Guo et al., 2004b) between two insoluble modes: accumulation and coarse (Fig. 1). The accumulation mode comprises 1.5 Mt of ash, and the coarse mode comprises 73.5 Mt of ash. Despite its small mass, ash in the accumulation mode is important since it has a much longer lifetime than ash in the coarse mode. We use the standard EMAC IMPORT and OFFEMIS submodels to initialize the SO₂, water vapor, and ash tracers (Kerkweg et al., 2006b). Section 1 in the Supplement explains the implementation of SO₂, water vapor, and volcanic ash injection mechanisms in EMAC.

In the main set of experiments, we release volcanic products in the specified model grid box centered on the altitude of 17, 20, or 25 km at the geographical coordinates of Mt. Pinatubo (15.1429 ^∘N, 120.3496 ^∘E) with the precalculated emission rates (in $\mathrm{molecules}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$) during 24 h. We call it a one-grid-box emission scheme (see Table 3 for details). In the 1s1-17km, 1s1-20km, and 1s1-25km experiments, we assume that only SO₂ is injected at 17, 20, and 25 km, respectively. In the 1w1-20km experiment we release SO₂ and water vapor at 20 km. The va0 experiments employ the same settings as 1w1 but assume injection of 75 Mt of ash. The va0 experiments do not account for the chemical aging of ash. The va1 experiments are similar to va0 but account for ash aging. In experiment 3s10-25km we inject SO₂ in the 3000 km wide latitude belt centered on the latitude of the eruption mimicking the setting in Brühl et al. (2015). The injected layer is 10 boxes thick (from 22.5 to 27.5km) and is centered on the altitude of 25 km. When this is not ambiguous, we refer to the clusters of experiments with the same physics using a generic name without specifying injection altitudes, such as 1s1, 1w1, va0, and va1. Experiments 1s1 are used to study the sensitivity to the height of the injection of volcanic SO₂. The 1w1 experiments with 150 Mt and 15 Mt injected water allows us to quantify the sensitivity to the mass of injected water vapor. Experiments va0 and va1 are designed to quantify the effect of ash and ash aging. Experiment va1-20km-12Mt is designed to study the sensitivity to the amount of injected SO₂.

To reduce the effect of internal model variability in each experiment, we calculate five ensemble members using different atmospheric initial conditions. The analysis in this study is performed and presented for the ensemble means. We show the “ensemble” variability for the selected variables.

5 Results

First, we compare the model results with observations, focusing on spatial–temporal distributions of SO₂, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, and other related chemicals. We also test the simulated stratospheric AOD (SAOD), which defines volcanic radiative effect, and the stratospheric temperature response, which measures volcanic climate impact against available observations. In addition, we compare the surface area density (SAD) that controls heterogeneous chemistry within the volcanic cloud and aerosol effective radius (R_eff), which characterizes aerosol size distribution (see Figs. 2–5). The spatially averaged R_eff is calculated as a ratio of the third M3_m and the second M2_m moments of each aerosol mode m integrated over the entire domain (Eqs. 1 and 2). The effective radii for individual modes and for the entire aerosol size distribution are given by Eqs. (3) and (4), respectively.

$$\begin{array}{}\text{(1)}& {\displaystyle }M{\mathrm{2}}_m=\underset{v}{\iiint }{N}_m{R}_{\mathrm{m}}^{\mathrm{2}}{\exp }^{\left(\mathrm{2}{\ln }^{\mathrm{2}}{\mathit{\sigma }}_m\right)}\mathrm{d}x\mathrm{d}y\mathrm{d}z,\text{(2)}& {\displaystyle }M{\mathrm{3}}_m=\underset{v}{\iiint }{N}_m{R}_{\mathrm{m}}^{\mathrm{3}}{\exp }^{\left(\frac{\mathrm{9}}{\mathrm{2}}{\ln }^{\mathrm{2}}{\mathit{\sigma }}_m\right)}\mathrm{d}x\mathrm{d}y\mathrm{d}z,\text{(3)}& {\displaystyle }{R}_{\text{eff}}^m={\displaystyle \frac{M{\mathrm{3}}_m}{M{\mathrm{2}}_m}},\text{(4)}& {\displaystyle }{R}_{\text{eff}}={\displaystyle \frac{{\sum }_{m=\mathrm{1}}^{m=N_{\text{modes}}}M{\mathrm{3}}_m}{{\sum }_{m=\mathrm{1}}^{m=N_{\text{modes}}}M{\mathrm{2}}_m}},\end{array}$$

where N_m is the number density for aerosol mode m, R_m is the median radius, and σ_m is the width of the aerosol mode m. N_modes is the number of aerosol modes.

Figure 2 (a) The 20^∘ S–20^∘ N equatorial average visible SAOD (0.55 µm) from the 1s1 experiment with the 17, 20, and 25 km emission heights, AVHRR 0.630 µm, scaled visible SAGE/ASAP, and 0.525 µm CMIP6. (b) SO₂ (solid lines) and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (dashed lines) globally integrated masses calculated using output from the 1s1 experiment with 17, 20, and 25 km emission heights; the observed Guo et al. (2004b) SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ masses are shown by markers. (c) Equatorial average effective radius from 1s1 experiment with 17, 20, and 25 km emission height, and SAGE/ASAP (“Stratospheric processes and their role in climate” published in the Assessment of Stratospheric Aerosol Properties) retrievals (solid black). Error bars represent 1σ from the ensemble mean.

Figure 3 Zonally average visible SAOD (0.55 µm) as a function of latitude and time: (a) 1s1-17km, (b) 1s1-20km, (c) 1s1-25km, (d) 3s10-25km, (e) scaled visible SAGE/ASAP, and (f) 0.630 µm AVHRR (Advanced Very High-Resolution Radiometer).

Figure 4 Zonally average surface area density (SAD, $\mathrm{\mu }{\mathrm{m}}^{-\mathrm{2}}{\mathrm{cm}}^{-\mathrm{3}}$) as a function of latitude and time at 20 km (a, d ,g and j), 25 km (b, e, h and k), and 30 km (c, f, i and l). The experiments are as follows: (a–c) 1s1-25km, (d–f) 1s1-20km, (g–i) 1s1-17km, and observed (j–l) SAGE/ASAP.

Figures 2–4 compare various parameters in 3s10-25km and 1s1 experiments with the different injection heights with the AVHRR and SAGE/ASAP observations. The AVHRR zonal mean visible SAOD is largely consistent in spatial–temporal behavior with the scaled SAGE/ASAP SAOD (Fig. 3). The original SAGE/ASAP visible SAOD is almost half of the AVHRR SAOD because the SAGE II sensor was saturated during the first few weeks after the eruption. Therefore data at the initial stage of eruption are sparse. The AVHRR sensed the entire atmospheric column, including the troposphere, the effect of which could be estimated only approximately (Thomason, 1992; Russell et al., 1996; Kremser et al., 2016). The consistency between the scaled SAGE/ASAP and AVHRR visible SAODs lessens in the late fall of 1991 when scaled SAGE/ASAP SAOD begins overestimating AVHRR SAOD. Discrepancies between different data sets are discussed in Bingen et al. (2004). Despite sparse observations at the initial stage of volcanic cloud development, SAGE/ASAP is the only global satellite observation that recorded the vertical structure of the Pinatubo cloud. For example, Fig. 4 demonstrates aerosol SAD at different altitudes as reported by SAGE/ASAP and simulated in the model.

Below we study the sensitivity of volcanic cloud evolution to all the main factors: injection height, amount of injected water, injection of ash, and ash aging. We start from sensitivity to injection height using the simplest 1s1 experiments with SO₂ only injections. The cloud height is essential because it defines the wind field that drives the dispersion of volcanic debris. The O₃ and water vapor mixing ratios, which affect chemical reactions and aerosol microphysics within the plume, are also height dependent.

5.1 Sensitivity to injection height

Figure 2a–c compare the observed and simulated SAOD, SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ masses and R_eff in the 1s1 experiments with different injection heights. The altitudes where volcanic debris resides depend not only on the initial injection height but on upward stratospheric motion and lofting driven by radiative heating of volcanic debris (Stenchikov, 2021; Niemeier et al., 2009; Kinnison et al., 1994; Aquila et al., 2012). The latter process and the rate of chemical transformations within a volcanic cloud are sensitive to the initial concentrations of optically and chemically active materials within a fresh volcanic cloud, i.e., in terms of our simulation settings, to the volume that a cloud initially occupies.

Experiment 3s10-25km assumes a zonally uniform SO₂ release at 25 km altitude within a latitude belt centered on the latitude of the eruption (15.1429^∘ N). In experiment 1s1-25km, we release SO₂ centered on the same height as in the 3s10-25km experiment but within one model grid box at the geographic coordinates of the Pinatubo eruption (15.1429^∘ N, 120.3496^∘ E). This causes initially higher SO₂ concentrations compared to the 3s10-25km experiment. The eruptive materials are released with a constant mass emission rate and spread for more than 1000 km during the 24 h of emission. Despite SO₂ having been released at the same altitude, these two experiments exhibit remarkable differences in the globally averaged SAODs (see Fig. 2a), ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ masses (Fig. 2b), and spatial distributions of SAOD and SAD (Figs. 3 and 4). To understand the mechanism of the strong sensitivity of the volcanic cloud evolution to its initial stage, we test the 3s10-25km experiment and the one-grid-cell 1s1 experiments with the 17, 20, and 25 km injection heights against observations.

5.1.1 SAOD

Contrary to the 1s1-25km experiment, the visible tropical SAOD in experiment 1s1-20km compares well with SAOD from the scaled SAGE/ASAP and AVHRR observations (Fig. 2a). The visible SAOD from the Sixth Coupled Model Intercomparison Project (CMIP6) (Eyring et al., 2016; Zanchettin et al., 2016) mimics the original visible SAGE/ASAP SAOD, develops slowly, and is half of the scaled SAGE/ASAP and AVHRR.

The equatorial average (20^∘ S–20^∘ N) SAOD in 1s1-17km is half the size of the 1s1-20km and 3s10-25km SAODs. The 1s1-25km SAOD is even smaller (Fig. 2a). All SAODs except that in the 1s1-25km experiment are bigger than the CMIP6 SAOD. The SAOD in the 3s10-25m experiment grows faster than the SAOD in the 1s1 experiments, reaching 0.33 in August 1991. At a given chemical composition of sulfate aerosol particles, the transient SAOD depends both on the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass, i.e., the rate of oxidation of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, and on aerosol size distribution, i.e., R_eff. The smaller sulfate aerosol particles have a bigger collective cross-section per unit mass than larger ones. So a bigger mass of large sulfate particles might have a smaller SAOD than a smaller mass of smaller sulfate particles. This must be considered when evaluating the mass of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ and the sulfate aerosol SAOD in observations and model experiments.

5.1.2 Oxidation of SO₂

Figure 2b shows the globally integrated SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ masses in the 3s10-25km and 1s1 experiments with the different emission heights as functions of time. The SO₂ mass in the 1s1-20km experiment decreases more slowly than in all other experiments. The ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass in the 3s10-25km grows faster than in the other experiments in Fig. 2b. This is because the SO₂ oxidation rate depends on the abundance of OH radicals. The OH production depends on O₃ concentration and incoming UV radiation. Because in the 3s10-25km experiment, SO₂ is distributed zonally over the entire latitude belt, its concentration in a volcanic cloud is lower than in all one-grid-box injection experiments. Hence, the SO₂ oxidation is more efficient in the 3s10-25km experiment than in the 1s1 runs because initially, there are more OH radicals available per SO₂ molecule in the latitude belt than in the one grid box. All 1s1 experiments underestimate ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass in the first few days in comparison with the available observations (Fig. 2b). The presence of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in a fresh volcanic plume detected in observations is puzzling as the models do not account for the physical mechanisms that could produce it in such a short time. To explain this discrepancy, Guo et al. (2004a) suggested that 1–2 Mt of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ was injected at the initial stage of the eruption. However, we do not account for the initial ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ release in this study.

5.1.3 Spatial–temporal evolution of SAOD and SAD

The spatial–temporal patterns of visible SAOD in the 1s1-20km experiment compare well with AVHRR and scaled SAGE/ASAP observations (Fig. 3). However, the poleward aerosol transport in the model is too fast. This is a known deficiency of global models, which simulate too transparent subtropical barriers due to coarse spatial resolution (Oman et al., 2006). The 1s1-25km visible SAOD is smaller than the scaled SAGE/ASAP and AVHRR SAODs and exhibits a qualitatively incorrect evolution of the volcanic cloud, which moves too far north, similar to that reported by Stenchikov et al. (2021) for volcanic injection at 24 km altitude. The 3s10-25km SAOD has a realistic spatial–temporal structure but substantially overestimates observed SAODs. SAOD in the 1s1-17km experiment (Fig. 3) exhibits faster poleward transport than in the 1s1-20km run due to more vigorous wave activity at lower altitudes in the stratosphere. In this experiment, the equatorial aerosol reservoir dissipates too quickly because of its proximity to the tropopause and intensive poleward transport.

Figure 4 compares the SAD in the 1s1-17km, 1s1-20km, and 1s1-25km experiments with the SAGE/ASAP observations (Thomason et al., 1997). SAD facilitates heterogeneous reactions in the volcanic cloud. Both sulfate aerosols and volcanic ash contribute to SAD, but in 1s1 experiments, we only account for sulfate aerosol surfaces. Therefore, the simulated SAD is expected to be smaller than the observed one, especially at the very beginning after the eruption. Only the 1s1-20km experiment shows SAD distributions consistent with observations at all three levels: 20, 25, and 30 km. In the 1s1-17km experiment the model underestimates SAD at 25 and 30 km. In the 1s1-25km experiment, the volcanic cloud climbs unrealistically high to 30 km and above. At that height, sulfate droplets tend to evaporate. The sulfuric acid photolyzes back to SO₂, which reaches the mesosphere (Rinsland et al., 1995).

5.1.4 Aerosol size distribution

Figure 2c compares R_eff from SAGE/ASAP averaged over the tropical belt, the 3S10-25km, and the 1s1 experiments with the 17, 20, and 25 km injection heights. In the control case, the model R_eff = 0.14 µm is lower than the observed unperturbed value of 0.17–0.19 µm (Russell et al., 1993), as the model underestimates the effect of anthropogenic sulfur emissions on the stratospheric Junge layer (Marandino et al., 2013; Brühl et al., 2015).

SAGE-II observations suggest that aerosol extinction increases, and its maximum shifts from 0.385 to 0.525 µm soon after the Pinatubo eruption, indicating the sudden increase in sizes of aerosol particles (Thomason, 1992; Thomason and Peter, 2006; Kremser et al., 2016). The observed effective radius increases from the background level to about 0.5 µm in 6 months (Russell et al., 1996).

In Fig. 2c, R_eff in the 1s1-20km experiment increases gradually, reaching maximum R_eff = 0.4 µm in September of 1991, and then decreases due to the settling of larger particles. In the 1s1-25km experiment, R_eff is the largest among all experiments. R_eff in this experiment initially grows faster than in all other runs, reaching maximum R_eff = 0.45 µm in August of 1991, and then decreases, merging with all other experiments in December of 1991. In the 1s1-17km experiment, R_eff is the smallest compared to other experiments, as it loses SO₂ mass through the tropopause. When testing our simulated R_eff against observations, we have to mention that the SAGE/ASAP R_eff is quite uncertain (Ansmann et al., 1997).

Figure 5 Perturbations (with respect to control) of 20^∘ S–20^∘ N equatorial average chemical constituents as a function of pressure (from 300 to 1 hPa) and time in 1s1-17km (a, d, g, j, m, p, and s), 1s1-20km (b, e, h, k, n, q, and t), and 1s1-25km (c, f, i, l, o, r, and u). (a–c) SO₂ (ppbv); (d–f) ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ (ng m⁻³); (g–i) OH (pptv); (j–l) H₂SO₄ (pptv); (m–o) NO_x (ppbv); (p–r) NO_y (ppbv), NO_x = (NO + NO₂), NO_y = (NO_x + NO₃ + HNO₃ + 2 × N₂O₅ + HONO + HNO₄ + ClONO₂ + BrONO₂); (s–u) O₃ (ppmv).

5.1.5 Impact on chemical composition

Figure 5 shows vertical cross-sections of the mixing ratios or concentrations of SO₂, ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$, OH, H₂SO₄, NO_x, NO_y, O₃ in the 1s1 experiments with the 17, 20, and 25 km injection heights averaged over the equatorial belt (20^∘ S–20^∘ N). We do not account for the SO₂ radiative effect (Stenchikov, 2021), and there is no ash injection in these experiments. Therefore, only sulfate aerosols cause heating and lofting of the volcanic cloud. Before a significant amount of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ develops, lifting is caused by regional upward motion in the Brewer–Dobson (B–D) circulation, which EMAC captures. The SO₂ and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ clouds separate due to gravitational settling of sulfate aerosols (Fig. 5a–f). This initiates multi-layer distributions of all other tracers. Stratospheric vertical uplift is getting stronger at higher altitudes (at least in EMAC). This is seen well in the 1s1-25km run in comparison with the 1s1-17km and 1s1-20km runs (see Fig. 5a–c).

The volcanic cloud in the 1s1-25km experiment rises to 30 km, significantly higher than in all other experiments. The SO₂ oxidation rate slows down as the temperature rises. Therefore, the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass is smaller than in the other experiments. In addition, in the 1s1-25km experiment, R_eff is higher than in other experiments. This factor tends to lower the SAOD since larger particles in the 1s1-25km experiment are less optically efficient per unit mass and have a lower lifetime with respect to gravitational settling. Therefore, SAOD and SAD in this experiment are smaller than in the other experiments in Figs. 2 and 4, respectively.

The experiments with different emission heights generate different SAOD (Fig. 2) and SAD (Fig. 4). This partially results from different SO₂ oxidation rates defined by the abundance of OH radicals at different altitudes. Three weeks after the eruption, OH is reduced around the injection height because of stratospheric water consumption by the oxidation of SO₂ in all three experiments in Fig. 5g–i. OH remains depleted above the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ cloud, where the SO₂ mixing ratio is high. The change in OH is generally the largest in the 1s1-25km experiment. The increase in H₂SO₄ is more pronounced in the 1s1-25km experiment (Fig. 5j–l). Initially, the H₂SO₄ increases at the emission level. This is seen until November of 1991. Then a secondary plume of H₂SO₄ is formed at a higher altitude, above 29 km. Oxidation of volcanic SO₂ in the stratosphere also perturbs the Chapman cycle and reduces the ozone mixing ratio in the stratosphere (Fig. 5s–u).

We account for 12 heterogeneous reactions. Heterogeneous reactions on aerosols explain the re-partitioning between NO_x and the reactive nitrogen reservoir NO_y. Following Danilin et al. (1999), we evaluate the effect of heterogeneous chemistry by the abundance of NO_x (NO + NO₂) and total inorganic nitrogen, NO_y (NO_x + NO₃ + HNO₃ + 2 × N₂O₅ + HONO + HNO₄ + ClONO₂ + BrONO₂). In addition, the heterogeneous chemistry might be influenced by the halogen reactions by chlorine and bromine activation. However, we do not account for the volcanic halogen injections in this study; therefore, their impact might be minor. Here, the main pathway in this transformation is the oxidation of NO_x to form N₂O₅, which interacts with water vapor to form HNO₃. N₂O₅ can also interact with halogens on the surface of aerosols (sulfate or ash). The heterogeneous reaction of N₂O₅ and water on the surface of aerosols effectively depletes NO₂ from the active reaction system depending on SAD (Fig. 4). During the daytime, HNO₃ can photo-dissolve and form OH, interacting with NO_x to form N₂O₅ and then HNO₃, while at nighttime, the formation of HNO₃ is one-way via oxidation of NO_x and N₂O₅. N₂O₅ can decompose back to NO₃ and NO₂ either photochemically or thermally, depending on the overhead column of ozone. As altitude increases, temperature increases, and the production rate of HNO₃ increases. The limiting factors in the heterogeneous formation of HNO₃ are that of NO₂, O₃ column, and SAD (Seinfeld and Pandis, 2006).

Figure 5m–r show the strong dependence of NO_x and NO_y on injection height within the aerosol cloud. The NO_x mixing ratio decreases, and the NO_y mixing ratio increases, along with the increase in the injection height. The changes in NO_x and NO_y affect the ozone cycle (Seinfeld and Pandis, 2006). The dependence of the background ozone concentrations on altitude adds to the sensitivity of the cloud evolution to injection height. Furthermore, the modified ozone concentration (Fig. 5s–u) feeds back to the OH production and hence the SO₂ oxidation rate. For the 1s1-17km experiment, the depletion of NO_x (Fig. 5m) is lower than for the 25 km injection (Fig. 5 o), while the production of NO_y at the 25 km (Fig. 5r) injection is higher than for the 17 km experiment (Fig. 5 p). At a higher altitude, the ozone concentration and SAD (Fig. 4c) is higher, and hence the formation of HNO₃ is higher for the 1s1-25km experiment (see Fig. S9 in the Supplement). Although the change in NO_y for 1s1 experiments at 17, 20, and 25 km injection is insignificant (see Fig. S9b in the Supplement), the heterogeneous transformation from N₂O₅ to HNO₃ is efficient. The transformation is enhanced (Fig. S9b–d in the Supplement) by the injection of ash particles due to the additional SAD and heating by ash and the associated stronger lofting of the volcanic plume.

Along with the chemical processes, the interaction of volcanic debris with the tropopause and the stratopause adds to the sensitivity of the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass to the height of the injection. In the 17 km injection height experiment, the cloud loses part of the mass through the tropopause, but in the 25 km injection height experiment, part of the sulfur is transported to the mesosphere and gets lost for immediate sulfate formation. It descends to the stratosphere again in high latitudes in winter. The volcanic debris injected at 20 km stabilizes in the middle of the stratosphere. Hence, in this case, the volcanic cloud is less affected by interaction with the tropopause and stratopause.

5.2 Water vapor intrusion

As expected, warming the tropical tropopause layer by radiative heating of volcanic debris facilitates the cross-tropopause troposphere-to-stratosphere transport of water vapor (Oltmans and Hofmann, 1995; Nedoluha et al., 1998; Joshi and Jones, 2009). The presence of extra water vapor in the stratosphere intensifies OH production and accelerates SO₂ oxidation (LeGrande et al., 2016).

Figure 6 Perturbations (with respect to control) of 20^∘ S–20^∘ N stratospheric (above 100 hPa) integrated masses (Mt) in the 1s1 experiments with injection heights 17, 20, and 25 km, as a function of time. (a) Water vapor; (b) ice; (c) aerosol water. Error bars represent 1σ from the ensemble mean.

For the 1s1-17km experiment, the stratospheric (i.e., above 100 hPa) water vapor mass integrated over the tropical belt increases by about 30 Mt (Fig. 6a). However, an increase in water vapor just above the tropopause does not affect volcanic cloud evolution much because the bulk of this water vapor is well below the altitude where the core of the volcanic cloud resides. Cross-tropopause water transport decreases as injection height increases. For example, the 1s1-25km experiment shows no cross-tropopause water transport. Three weeks after the injection, the aerosol water associated with sulfate aerosols in the 1s1-17km experiment (Fig. 6c) is higher than in other experiments because, in the 1s1-17km experiment, the volcanic cloud is closer to the tropopause than in the other experiments. In the 1s1-20km run, the mass of ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ continues increasing during August and September 1991 (Fig. 2b), and the associated aerosol water also increases to 3.5 Mt as shown in Fig. 6c. Little ice is accumulated in the stratosphere in all experiments (Fig. 6b) since it is quickly removed by gravitational sedimentation. In the 1s1-20km experiment liquid water mass peaks at 3 Mt, in the 1s1-17km at 2 Mt, and in 1s1-25km at 1 Mt (Fig. 6d).

5.3 Volcanic water injection

Water vapor injected into the stratosphere with a volcanic plume could directly affect the initial evolution of a volcanic cloud since it is concentrated within it. Most of this water is brought by the entrainment of tropospheric water in an explosive jet or co-ignimbrite convective updrafts; nevertheless, the term volcanic water is used here. A wide range (75–150 Mt) of volcanic water vapor injection for the Pinatubo eruption was reported (Joshi and Jones, 2009; Nedoluha et al., 1998). However, the amount of volcanic water retained in the stratosphere depends on the injection height. Almost all water vapor injected at a low temperature just above the tropopause forms ice and quickly sediments (Stenchikov et al., 2021). A larger fraction of water vapor injected at higher altitudes, with higher stratospheric temperatures, could remain in the stratosphere. To test the sensitivity of volcanic clouds to the amount of volcanic water, we conduct the 1w1 simulations injecting simultaneously SO₂ and different amounts of water vapor.

Figure 7 Perturbations (with respect to control) of 20^∘ S–20^∘ N stratospheric (above 100 hPa) SAOD and integrated masses (Mt) in the 1w1-20km and 1w1-25km experiments with respect to, correspondingly, 1s1-20km and 1s1-25km experiments as a function of time. Left column (a, c, e, g and i): 1w1 experiments with 15 Mt of volcanic water vapor injection. Right column (b, d, f, h and j): 1w1 experiments with 150 Mt water vapor injection. (a, b) SAOD; (c, d) ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$; (e, f) aerosol water; (g, h) ice; (i, j) water vapor.

Figure 7 compares the time series of the equatorial SAODs and changes in the globally integrated masses of sulfate and water species in the stratosphere (above 100 hPa) in the 1w1 experiments, with the simultaneous injection of SO₂ and 15 or 150 Mt of volcanic water vapor at 20 and 25 km with the corresponding 1s1 experiments (see Table 3). Water species comprise water vapor, ice, and aerosol water. The aerosol water accumulates in sulfate and ash particles.

The effect of volcanic water on the generation of the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass and SAOD depends on the amount of water vapor retained in the stratosphere after the eruption. The sensitivity of SAOD and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass to the injected volcanic water vapor is higher in the 1w1-25km experiment compared to the 1s1-20km experiment (Fig. 7c and d). The increase in sulfate mass results from the acceleration of SO₂ oxidation facilitated by the higher water vapor concentration (Seinfeld and Pandis, 2006). The sensitivity to the amount of water vapor in the experiments with the 20 km injection altitude is weaker than in the experiments with the 25 km injection height because most of the water vapor injected at 20 km condenses and deposits from the stratosphere since the temperature at 20 km is lower than at 25 km (Fig. 7i).

Figure 8 Globally integrated stratospheric volcanic ash mass (Mt) above 70 hPa as a function of time in the va0-20km and va1-20km experiments. (a) Dry accumulation mode. (b) Dry coarse mode. (c) Wet accumulation mode. (d) Wet coarse mode. Error bars represent 1σ from the ensemble mean.

Figure 9 Globally integrated stratospheric masses (Mt) as a function of time. (a) Volcanic ash in va0-20km, va1-20km as well as in AVHRR, and HIRS (High-Resolution Infrared Radiation Sounder) retrievals (Guo et al., 2004b). (b) SO₂ in the 1s1-20km, 1w1-20km, va0-20km, and va1-20km experiments as well as in TOMS (Total Ozone Mapping Spectrometer) and TOVS (TIROS Operational Vertical Sounder) observations.

5.4 Volcanic ash injection

In the va0 and va1 experiments, we inject 75 Mt of ash together with SO₂ and water vapor. The va1 experiment accounts for ash aging, but the va0 does not. In both experiments, we assume that volcanic ash is initially hydrophobic. Therefore, we inject it into the insoluble (dry) accumulation and coarse modes (Figs. 1 and 8). In the va1 experiments, volcanic ash ages quickly, populating the soluble (wet) modes (Fig. 8c and d), while ash particles in the va0 experiments remain in the dry modes (accumulation and coarse). In the va1 experiments, ash particles increase in size due to the aging and associated water and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ uptake, which tends to transfer particles from the accumulation to the coarse mode.

In the va0 experiments, ash in the coarse mode (see Fig. 8b) sediments from the stratosphere for 2 d, but ash particles in the accumulation mode remain in the stratosphere for 1 week (Fig. 8a). In the va1 experiments, the ash mass in the wet modes increases quickly due to dry-to-wet particle conversion shown by the arrow in Fig. 8. The aging of ash particles slows the decrease in ash mass in both accumulation and coarse modes.

In the experiments with a 25 km injection height, it takes longer for ash to reach the tropopause and leave the stratosphere than in the experiments with the 20 km injection height. For instance, after the first day of injection, 60 Mt of insoluble coarse ash remains in the stratosphere (not shown) for the va0 experiment with 25 km injection height compared to 1.7 Mt when ash is injected at 20 km (see Fig. 8b).

Figure 10 Visible SAOD (0.55 µm), aerosol effective radius, R_eff above 100 hPa, and integrated masses (Mt) simulated in the 1s1-20km, 1w1-20km, va0-20km, and va1-20km experiments as a function of time. (a) Simulated, as well as observed AVHRR, scaled SAGE/ASAP, and CMIP6 20^∘ S–20^∘ N SAODs. (b) Simulated stratospheric (above 100 hPa) R_eff as well as that observed in SAGE/ASAP in 20^∘ S–20^∘ N. (c) Stratospheric ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in accumulation mode in 20^∘ S–20^∘ N. (d) Stratospheric ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in coarse mode in 20^∘ S–20^∘ N. (e) Stratospheric aerosol water in accumulation mode in 20^∘ S–20^∘ N. (f) Stratospheric aerosol water in coarse mode in 20^∘ S–20^∘ N. (e) Tropospheric (below 100 hPa) ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in accumulation mode integrated globally. (f) Tropospheric ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in coarse mode integrated globally.

Figure 11 The 20^∘ S–20^∘ N average perturbations of chemical constituents as a function of pressure and time in va0-20km and va1-20km experiments. (a) SO₂ in va0-ctr (ppbv). (b) SO₂ in va1-ctr (ppbv). (c) ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in va0-ctr (ng m⁻³). (d) ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ in va1-ctr (ng m⁻³). (e) H₂SO₄ in va0-1w1 (ppbv). (f) H₂SO₄ in va1-va0 (ppbv). (g) OH in va0-1w1 (pptv). (h) OH in va1-va0 (pptv). The contour lines show the accumulation mode ash mixing ratio (ppbv); orange contour lines for va0 and blue contour lines for va1-va0 in (e–h).

Figure 12 Shortwave heating rate (K d⁻¹) averaged over the tropical belt (20^∘ S–20^∘ N) shown as a function of time and height overlaid by contours of volcanic ash mass concentration (accumulation mode) for (a) va0-20km and (b) va1-20km experiments. Longwave heating rate (K d⁻¹) for (c) va0-20km, and (d) va1-20km experiments, both overlaid by contours of volcanic ash mass concentration (coarse mode). All heating rates are calculated by a double call of the radiation routines.

Figure 9a shows the evolution of the stratospheric ash masses in the va0 and va1 experiments compared to the AVHRR and HIRS/2 retrievals (Guo et al., 2004b). In the va1 experiment, the model ash mass is higher than in the AVHRR and HIRS/2 observations, while in the va0 experiment, the ash mass is underestimated compared with observations. However, the uncertainties in the AVHRR-derived ash mass are ± 53 % (Gu et al., 2003) and ± 85 % in HIRS retrievals (Yu and Rose, 2000). The larger ash mass in the va1 experiment (compared to va0) on the first day after the eruption causes stronger heating and lofting of ash in the va1 experiment than in va0, which prolongs its lifetime in the stratosphere (Fig. 12).

The injection of volcanic ash significantly increases the stratospheric optical depth and R_eff during the few days after the eruption. This is shown in Fig. 10a and b, which compare the time series of SAOD and the effective radius averaged within 20^∘ S–20^∘ N latitude belt above 100 hPa in the 1s1-20km, 1w1-20km, va0-20km, and va1-20km experiments, with available observations. The AVHRR and scaled SAGE/ASAP SAODs are consistent for 4–5 months after the eruption. The CMIP6 SAOD appears to be half the size when compared to them. During the 4–5 months following the eruption, the simulated SAOD (Fig. 10a) is slightly larger than in observations but decreases more quickly than in observations later on in all experiments except va0-20km. The SAODs in the va1-20km and va0-20km experiments grow more rapidly during the first 2 months than in the 1s1-20km and 1w1-20km experiments (Fig. 10a). The effective radii in the va1 and va0 experiments spike to about 0.6–0.8 µm during the first week after the eruption, when a significant amount of ash is present in the volcanic cloud.

Figure 10c and d show the evolution of the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass in the coarse and accumulation modes integrated over the 20^∘ S–20^∘ N latitude belt. In the va0 and va1 experiments, the stratospheric sulfate mass increases more rapidly than in the 1w1 and 1s1 experiments. This is consistent with SAOD in Fig. 10a and with the more rapid depletion of SO₂ mass in Fig. 9b, which demonstrates a better agreement with SO₂ mass observations. We relate the faster SO₂ oxidation in the va0 and va1 experiments to the effect of heterogeneous reactions on ash particles and more intensive volcanic cloud dispersion facilitated by ash radiative heating.

Two months after the eruption, in the 1s1 and 1w1 experiments, the sulfate mass in accumulation and coarse modes reaches maximums of 9 and 0.7 Mt, respectively. Thus, the sulfate formation rate increases in the va1 and va0 experiments compared to experiments without ash in both accumulation and coarse modes. The ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass reaches the maximum 2 weeks earlier in experiments with ash than in experiments without ash (Fig. 10c and d).

The aerosol water mass increases when sulfate mass increases, both in the accumulation and the coarse modes (Fig. 10e and f). A sulfate mass of 9 Mt is associated with an aerosol water mass of 3 Mt in the accumulation mode (Fig. 10c and e). This is consistent with the 75 % sulfuric acid solution assumed by Stenchikov et al. (1998). For the coarse mode, the aerosol water mass of 0.5 Mt is associated with roughly 0.8 Mt of sulfate (Fig. 10d and f). Both sulfate and wet ash particles accumulate aerosol water. In the long run, due to the shorter lifetime of ash particles, aerosol water is mainly retained in sulfate aerosols. Figure 10g and h show ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass in the coarse and accumulation modes in the troposphere (integrated below 100 hPa) globally for the same experiments. The tropospheric ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass of volcanic origin comprises ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ sedimented from the stratosphere. Because of the rapid wet removal, little sulfate (not exceeding 0.4 Mt in each mode) is accumulated in the troposphere both in accumulation and coarse modes. This is more than 1 order of magnitude less than the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass in the stratosphere. More sulfate mass sediments into the troposphere in the va0 and va1 experiments than in the 1s1 and 1w1 runs (Fig. 10e) because in the va1 and va0 experiments, the stratospheric sulfate mass is larger than in the 1s1 and 1w1 experiments.

5.5 Ash aging

Ash particles provide surface areas, enhancing the heterogeneous reactions and leading to significant changes in stratospheric chemistry (Danilin et al., 1999). Ash SAD is especially important in the first week after the eruption when a few sulfate aerosols form. Figure 11a–d show the 20^∘ S–20^∘ N mean SO₂ mixing ratio and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ concentration as a function of time and height for the va0-20km and va1-20km experiments. Ash radiative heating causes lofting of the SO₂ plume by about 1 km d⁻¹ in both experiments, similar to that found in Stenchikov et al. (2021), although ash in our simulations is more absorbing than in Stenchikov et al. (2021).

In both cases, the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ cloud is below 35 km, but SO₂ reaches the stratopause. Therefore some SO₂ penetrates the mesosphere. This effect is more significant in the va0 experiment because of slower SO₂ oxidation compared with the va1 experiment. The enhanced mixing ratio of SO₂ in the mesosphere above 45 km was detected in ATMOS observations (Rinsland et al., 1995) and simulated by Brühl et al. (2015) after the 1991 Pinatubo eruption.

Figure 11e and g show the changes in mixing ratio of H₂SO₄, and OH concentrations in the va0-20km experiment with respect to 1w1-20km, demonstrating the effect of ash injection. Figure 11f and h show changes in the same characteristics, except in the va1-20km experiment with respect to the va0-20km experiment to demonstrate the effect of aging. If aging is turned on, H₂SO₄ condenses on volcanic ash, decreasing sulfuric acid concentration. At the same time, ash facilitates heterogeneous reactions. Moreover, the stratospheric water vapor mass is larger for the va1 experiment than va0 because of the more substantial heating (longwave and shortwave) in the va1 experiment (Fig. 11b and c). The more significant heating of the va1 experiment increases the tropospheric water vapor transport into the stratosphere. As a result, the stratospheric water vapor mass in va1 is 15 Mt higher than in va0 after 4 months of the eruption. The increased stratospheric water vapor facilitates OH formation (11h). The combined effect of ash radiative heating, aging, and heterogeneous chemistry in our setting resulted in an increase in sulfuric acid and sulfate mass by about 10 %–20 %, compared to those experiments without ash injections.

Figure 12 shows the shortwave and longwave heating rates averaged over the tropical belt (20^∘ S–20^∘ N) caused by volcanic cloud for the va0-20km (panels a and c) and va1-20km (panels b and d) experiments, as a function of time and height. The contour lines show the ash concentrations for the accumulation (panels a and b) and coarse (panels c and d) modes. Both experiments have two distinguished periods of increased radiative heating (Fig. 12a and b). The first time the heating peaks just after the eruption and the second time 10 d later. The first peak is associated with ash absorption and the second with sulfate aerosol absorption. In both cases the SW heating by sulfate peaks at 25 km (see Fig. 11a and b). The average ash SW heating rates is about 0.4 K d⁻¹ in experiment va1-20km and 0.15 K d⁻¹ – in experiment va0-20km (Fig. 12a and b). The IR radiation cools the top of the volcanic cloud during the first few days after the eruption. Still, the absorption of upward IR radiation heats the bottom of the volcanic cloud. Heating caused by the absorption of IR radiation by sulfate aerosols is seen in about 10 d when enough ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ is generated. Figure 12d shows that in the va1-20km experiment, the longwave heating rate reaches 0.2 K d⁻¹. Radiative heating in the va0-20km experiment is weaker than in the va1-20km experiment (Fig. 12c). Thus, aging significantly enhances the radiative effect of ash for about 1 week after the eruption.

5.6 Long-term climate response to volcanic forcing

In Sect. 5.1, we showed that the model demonstrates strong SAOD sensitivity to the injection height during the first 6 months after the eruption. We also found that the spatial–temporal SAOD development in simulations with the volcanic emission of 17 Mt SO₂ at 20 km best fits the observations 6 months after the eruption, but SAOD is overestimated. Here we further test the volcanic cloud evolution and stratospheric temperature response for the post-Pinatubo period of 2.5 years against observations. We take advantage of the fact that the climate response provides another constraint on SAOD (Stenchikov et al., 1998; Kirchner et al., 1999). We also quantify the sensitivity of volcanic cloud evolution to the amount of injected SO₂ considering the 12 Mt SO₂ emission at 20 km height.

Figure 13 Visible and near-IR SAODs in the va1-20km and va1-20km-12Mt experiments as well as in AVHRR, scaled SAGE/ASAP, and CMIP6. (a) Visible 20^∘ S–20^∘ N SAODs as a function of time, (b) near-IR 20^∘ S–20^∘ N SAODs as a function of time, (c) visible globally averaged SAODs as a function of time, (d) near-IR globally averaged SAODs as a function of time, (e) simulated visible zonally average SAOD in va1-20km as a function of time and latitude, (f) simulated near-IR zonally average SAOD in va1-20km as a function of time and latitude, (g) SAGE/ASAP scaled visible zonally average SAOD as a function of time and latitude, and (h) SAGE/ASAP near-IR zonally average SAOD as a function of time and latitude. Error bars represent 1σ from the ensemble mean.

Figure 13 compares the post-eruption evolution of SAODs (visible and near-IR) and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass in the va1 experiments with the 20 km injection height and the 17 Mt and 12 Mt SO₂ emission with the observations from CMIP6, AVHRR, and SAGE/ASAP (scaled visible and original near-IR) SAODs for 2.5 years. In the va1-20km experiment with the 17 Mt SO₂ injection, the simulated SAOD overestimated the AVHRR and scaled SAGE/ASAP SAODs both in the visible and near-IR range (Fig. 14a–d) in July–August 1991. In the experiment with the 12 Mt SO₂ emission, the SAOD reduces and fits the observations in visible and near-IR range. It is important that the initial rate of development of visible and near-IR SAODs is similar in the model and observations both in the tropics and globally. It suggests that, at the stage when the aerosol cloud is still confined to the tropics, the model captures the SO₂ oxidation process and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ development quite well.

Starting from September 1991, the excessive speed of the poleward transport of aerosols in the model causes a faster decrease in SAOD in the simulations (both in the tropics and globally) than in the observations (Fig. 13e–h). This is because sulfate aerosols are mainly deposited in the midlatitude storm tracks through tropopause faults and in the polar regions in the downward branch of B–D circulation (Gao et al., 2007). The faster poleward aerosol transport facilitates both these processes.

Figure 14 Atmospheric temperature anomalies (K) for the post-Pinatubo period with respect to the 1990–2000 climatology from the va1-20km (a, b), va1-20km-12Mt (c, d), and MERRA-2 reanalysis (e, f). Panels (a, c, e) depict zonally average anomalies at the 50 hPa pressure level as a function of time and latitude, and (b, d, f) depict globally (70^∘ S–70^∘ N) averaged anomalies as a function of time and height/pressure.

We further evaluate the long-term model stratospheric temperature response to test the consistency of va1 simulations with observations. Figure 14 shows the temperature anomalies in the 1s1, va1-20km, and va1-20km-12Mt experiments and MERRA2 reanalysis. Figure 14a, c, and e depict the Hovmöller diagrams of zonal mean temperature anomaly at 50 hPa, and Fig. 14b, d, and f are the global mean (70^∘ S–70^∘ N) temperature anomaly as a function of height (or pressure) and time. All experiments in Fig. 14 resemble the spatial–temporal structure of the stratospheric temperature response well. They reproduce stratospheric volcanic heating in the first year after the eruption and the additional heating associated with the change in the QBO phase in 1993.

In the va1-20km-12Mt experiment, the SAOD is about 30 % lower than in va1-20km, causing lower radiative heating. In the va1-20km experiment, the peak of temperature response reaches 4 K at 50 hPa. In the va1-20km-12Mt experiment, the temperature peak is about 2.5 K (Fig. 14a and c), which better agrees with observations (Labitzke and McCormick, 1992) and the MERRA2 reanalysis (Fig. 14e). Figure 14b, d, and f show a peak temperature anomaly at 30 hPa in the model simulations and the reanalysis. Again, the temperature response in the va1-20km-12Mt experiment fits the MERRA2 temperature anomalies better than the va1-20km experiment (Fig. 14b, d, and f). Thus, reducing the injected SO₂ mass to 12 Mt shows a better agreement with the observations (Labitzke and McCormick, 1992), and the MERRA2 temperature anomalies, resulting in more realistic heating at 50 hPa in both the tropics and the subtropics.

6 Conclusions

In this study, we use the EMAC model with well-developed stratospheric chemistry (including heterogeneous chemistry) and detailed aerosol microphysics to explore the evolution of the volcanic cloud from the 1991 Pinatubo eruption, the largest and best-observed volcanic event in the 20th century. We tested the model results with observations of volcanic clouds and their radiative effect. We conducted ensemble simulations to study the volcanic cloud sensitivity to the injection height and the initial volume where volcanic materials were released (one grid box versus a latitude belt), as well as considering the effects of the co-injection of water vapor, ash, and ash aging on the formation of the volcanic cloud.

The model simulations with 20 km injection height exhibit the best agreement of the spatial–temporal evolution of SAOD with the AVHRR SAOD and SAGE/ASAP SAOD and SAD. In the 20 km injection experiments, the volcanic cloud is lifted to a height of 25 km by radiative heating, while in the experiments with volcanic material injection at 25 km, the volcanic cloud overshoots the 30 km level. The vertical distribution of SAOD and SAD in the observations and the model experiments with the 20 km injection height shows that the aerosol cloud stabilizes in the middle of the stratosphere at 25 km. In the experiments with the 17 and 25 km injection heights, the volcanic cloud interacts with the tropopause and the stratopause, respectively, causing some aerosol mass to escape from the stratosphere. The stratospheric oxidation capacity and wind fields are different at different altitudes, strengthening the sensitivity to the injection height. In the experiments with the zonally uniform SO₂ injection in a latitude belt at the height of 25 km, lofting is weak, so the volcanic cloud remains at the same altitude. The SAOD in the tropical belt in this experiment is overestimated due to the higher oxidation rate.

Because of the coarse spatial resolution (T42L90), similar to other global models, EMAC simulates too fast aerosol poleward transport with too quick an escape of the volcanic materials from the tropical stratosphere, accelerating the loss of the aerosol mass to deposition at the poles, in tropopause folds, and storm tracks.

The increase in water vapor in the stratosphere leads to the rise in the oxidation rate of SO₂ to ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$. The water vapor could be brought into the stratosphere by an eruptive jet, co-ignimbrite convection, and/or intruded through the tropopause heated by the absorption of solar and IR radiation by volcanic debris. The cross-tropopause water vapor intrusion does not affect the volcanic cloud evolution much, as most of the water penetrating through the tropopause accumulates below the volcanic cloud. The water vapor directly injected into the volcanic cloud in the 1s1-20km experiment increases the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass and SAOD by about 5 %. The sensitivity of the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass to the amount of injected water in the experiment with the 20 km injection height is low because most of the water vapor freezes and is quickly removed from the stratosphere in agreement with Stenchikov et al. (2021). So the masses of the remaining stratospheric water vapor in the 1w1-20km experiments with 15 and 150 Mt of water vapor injections do not differ much. A significant acceleration of SO₂ oxidation due to the injection of water vapor (LeGrande et al., 2016) is only reproduced in the experiments with the 25 km injection height when a significant mass of injected water vapor is retained in the stratosphere.

Our experiments' simulated mass of ash is within the AVHRR and HIRS estimates, but observations themselves are uncertain. Volcanic ash provides SAD for heterogeneous chemistry. This is most important during the first few weeks after the eruption when ash is still abundant but sulfate aerosol is not yet developed. The simultaneous injection of water vapor and non-aging ash in the va0-20km experiment increases the maximum SAOD and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass by 10 %.

In the va1-20km experiment, ash particles in the accumulation and coarse modes are entirely aged within 1 d after the injection. Aging increases the mass of ash particles. They continue uptaking water and ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ molecules until removed by transport or sedimentation. The coarse ash particles deposit within 1 week, while it takes 6 months to reduce the mass of the ash accumulation mode from 1.2 to 0.3 Mt. Overall, aging increases the SAOD by 20 % and the ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$ mass by 10 %. Aging increases the radiative effect of ash both in SW and IR. The injections of volcanic water vapor and ash significantly accelerate the formation of the sulfate aerosols during the first 2 months after the eruption in the va1-20km and va0-20km experiments.

The simulated maximum SAOD and stratospheric temperature anomalies in the va1-20km-12Mt experiment with the 12 Mt SO₂ injection closely resemble the temperature anomalies obtained from the reanalysis both in latitude and height. The inclusion of volcanic ash adds to the radiative heating of the volcanic debris during the first week after the eruption in agreement with Stenchikov et al. (2021), showing that the initial local heating results in the lofting of the aerosol cloud. Our simulations show that the interactive calculations of OH and heterogeneous chemistry increase the volcanic cloud sensitivity to water vapor and ash injections and have to be accounted for in simulations of volcanic impacts on climate and stratospheric chemistry.