Interactive stratospheric aerosol models' response to different amounts and altitudes of SO₂ injection during the 1991 Pinatubo eruption

A previous model intercomparison of the Tambora aerosol cloud has highlighted substantial differences among simulated volcanic aerosol properties in the pre-industrial stratosphere and has led to questions about the applicability of global aerosol models for large-magnitude explosive eruptions prior to the observational period. Here, we compare the evolution of the stratospheric aerosol cloud following the well-observed June 1991 Mt. Pinatubo eruption simulated with six interactive stratospheric aerosol microphysics models to a range of observational data sets.

Our primary focus is on the uncertainties regarding initial SO₂ emission following the Pinatubo eruption, as prescribed in the Historical Eruptions SO₂ Emission Assessment experiments (HErSEA), in the framework of the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA-MIP). Six global models with interactive aerosol microphysics took part in this study: ECHAM6-SALSA, EMAC, ECHAM5-HAM, SOCOL-AERv2, ULAQ-CCM, and UM-UKCA. Model simulations are performed by varying the SO₂ injection amount (ranging between 5 and 10 Tg S) and the altitude of injection (between 18–25 km).

The comparisons show that all models consistently demonstrate faster reduction from the peak in sulfate mass burden in the tropical stratosphere. Most models also show a stronger transport towards the extratropics in the Northern Hemisphere, at the expense of the observed tropical confinement, suggesting a much weaker subtropical barrier in all the models, which results in a shorter e-folding time compared to the observations. Furthermore, simulations in which more than 5 Tg S in the form of SO₂ is injected show an initial overestimation of the sulfate burden in the tropics and, in some models, in the Northern Hemisphere and a large surface area density a few months after the eruption compared to the values measured in the tropics and the in situ measurements over Laramie. This draws attention to the importance of including processes such as the ash injection for the removal of the initial SO₂ and aerosol lofting through local heating.

1 Introduction

Large-magnitude volcanic eruptions can emit sulfur dioxide (SO₂) and other gases directly into the stratosphere. An abrupt increase in stratospheric SO₂ creates a long-lived volcanic aerosol cloud that scatters incoming solar radiation, absorbs solar and infrared radiation, and affects the composition of the stratosphere. Such volcanically induced enhancements of the stratospheric aerosol layer exert strong direct effects on climate because they influence the Earth radiation budget and cool the surface via the reduced insolation (McCormick et al., 1995; Soden et al., 2002); they also show a range of indirect effects, due to the volcanic aerosols effects on stratospheric circulation, dynamics, and chemistry (e.g., Robock et al., 2009; Timmreck et al., 2012; Kremser et al., 2016).

Here we investigate the evolution of the volcanic aerosol cloud after the Mt. Pinatubo eruption in June 1991 by analyzing coordinated simulations within the HErSEA (Historical Eruptions SO₂ Emission Assessment) experiments, in the framework of the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA-MIP; Timmreck et al., 2018). Mt. Pinatubo is located in the western part of the island of Luzon, Philippines (15.1^∘ N, 120.4^∘ E). After preliminary eruptions from 12 June 1991, the climatic phase started at 05:30 UTC on 15 June 1991 and lasted for approximately 9 h. The volcanic cloud contained gases and particles of ice, ash, and sulfate and reached a maximum altitude of 40 km (Holasek et al., 1996). Ice and ash burden peaked at about 80 and 50 Tg, respectively, and early-formed sulfate mass was estimated at 4 Tg, based on infrared satellite data from the Advanced Very High Resolution Radiometer (AVHRR), the TIROS Operational Vertical Sounder (TOVS), and High Resolution Infrared Radiation Sounder/2 (HIRS/2) sensors (Guo et al., 2004a). Initial sulfur dioxide (SO₂) mass estimates from the ultraviolet Total Ozone Mapping Spectrometer (TOMS) and infrared TOVS sensors, indicated that the eruption injected 14–22 Tg of SO₂ (Bluth et al., 1992; Guo et al., 2004a). Other uncertainties pertain to the vertical extension of the volcanic cloud: SO₂ mass was injected between 18–30 km (Bluth et al., 1992; Baran et al., 1993) and concentrated around 25 km, over a rich ash layer peaking around 22 km (Guo et al., 2004b). The sulfate aerosol cloud peaked at 14 Tg in September (Lambert et al., 1993; Baran and Foot, 1994), with the largest aerosol concentration between 20 and 25 km of altitude and much lower amounts between 15 and 20 km (Winker and Osborn, 1992a, b; DeFoor et al., 1992). Recent volcanic SO₂ emission databases suggest for Pinatubo an amount and location of SO₂ emitted between 15 and 18 Tg of SO₂, at an altitude of between 19 and 28 km (Independent Volcanic Eruption Source Parameter Archive Version 1.0, ivespa.co.uk, VolcanEESM: Global volcanic sulphur dioxide (SO₂) emissions database from 1850 to present - Version 1.0, Multi-Decadal Sulfur Dioxide Climatology from Satellite Instruments; Aubry et al., 2021; Neely III and Schmidt, 2016; Carn, 2022).

Several modeling studies have evaluated the simulated global and tropical sulfate loadings compared to observations, with some studies (Niemeier et al., 2009; Toohey et al., 2011; Brühl et al., 2015) finding agreement when emitting in the mid-range of the best-estimate stratospheric SO₂ loading of 14–22 Tg SO₂ (Guo et al., 2004a). In contrast, a number of recent studies found agreement only when injecting an amount of SO₂ below the lower limit observed of 10 Tg SO₂, considering different injection heights and vertical distributions (Dhomse et al., 2014; Sheng et al., 2015a; Mills et al., 2016); this difference partly motivates the design of the ISA-MIP HErSEA intercomparison (see Timmreck et al., 2018). Approaching the problem from a model intercomparison perspective, different past projects have revealed large differences in the simulation of the aerosol radiative forcing, and not just for Pinatubo.

A first multi-model intercomparison study of global stratospheric interactive aerosol models was set up in the frame of the Model Intercomparison Project on the climatic response to Volcanic forcing (VolMIP; Zanchettin et al., 2016). To create a common forcing data set for VolMIP experiments which considers a volcanic eruption with radiative forcing comparable to that of the 1815 Tambora eruption, a pre-study was set up (Marshall et al., 2018). This VolMIP-Tambora ISA experiment establishes a well-defined set of injection parameters to simulate the Tambora volcanic aerosol cloud interactively with stratospheric aerosol models. Multi-model analysis of the simulated volcanic aerosol distribution shows large inter-model differences (Marshall et al., 2018; Clyne et al., 2021).

Marshall et al. (2018) used Arctic and Antarctic ice core information about sulfate deposition to constrain the VolMIP-Tambora ISA model simulations. The four models involved in this experiment revealed large discrepancies in the simulated aerosol burden (50–58 Tg SO₄ at the peak), resulting in deposition magnitudes in Antarctica ranging from 19 to 264 kg km⁻². They attributed the differences between the models, and between models and observations, to different sulfate formation and transport through meridional circulation and stratosphere–troposphere exchange and different deposition schemes. The contribution to the overall uncertainty of the sulfate formation processes was then further investigated in a subsequent study by Clyne et al. (2021), which focused on the evolution of the global stratospheric aerosol optical depth. The reasons for the discrepancies between the models were attributed to differences in particle size, which influence the scattering efficiency and the lifetime of the stratospheric aerosols and the treatment of hydroxyl radical (OH) chemistry, which in turn affects the timing of sulfate formation.

The Geoengineering Model Intercomparison Project Phase 6 (GeoMIP6; Kravitz et al., 2015) also includes experiments with injection of stratospheric sulfate aerosol precursors (G6Sulfur) in an amount necessary to reduce the net radiative forcing from the SSP5-8.5 scenario to the SSP2-4.5 one. Participating models in G6Sulfur directly injected SO₂ in the tropical stratosphere with different altitude and latitude ranges of injection or prescribed the aerosol optical depth or aerosol distribution derived from previous simulations. The amount of SO₂ required to achieve the proposed cooling varies by a factor of 2 between models and results in a different temporal and latitudinal distribution of aerosols that affects surface temperature and local precipitation differently (Visioni et al., 2021).

In contrast to the aforementioned model intercomparison studies, the ISA-MIP HErSEA experiments offer a test of the reliability of these models by allowing a direct comparison of the simulated volcanic enhancement of the stratospheric aerosol layer with observation data sets, especially during the Mt. Pinatubo eruption, for which several satellite and in situ measurements are available. Hence, HErSEA was developed to determine which set of volcanic emission source parameters allows models to reproduce the available measurements, and understand how their different chemical and microphysical schemes, stratospheric dynamics, and radiative transfer treatment influence these choices. Specifically, HErSEA focuses on the uncertainty in the initial volcanic emission in terms of amount and injection altitude of SO₂ for the recent large-magnitude volcanic eruptions in the last 100 years (Mt. Agung in 1963, Mt. El Chichón in 1982, Mt. Pinatubo in 1991); multiple interactive stratospheric aerosol simulations of each of the volcanic aerosol clouds with common upper-, mid-, and lower-estimate amounts and injection altitudes of sulfur dioxide were performed. Here we investigate the evolution of the volcanic aerosol cloud after the Mt. Pinatubo eruption by analyzing Atmospheric Model Intercomparison Project (AMIP)-type (Gates et al., 1999) simulations within the HErSEA framework. In particular, we ask whether previous results in inter-model differences are confirmed in this new MIP; the presence of multiple injection settings common between all models will also allow an exploration of the reason for these differences, based on the models' abilities to reproduce observations with different sets of initial conditions of the volcanic emissions.

The experimental design, the main features of the participating models, and the observational data sets are described in Sect. 2. Section 3 shows model results of the optical and microphysical properties of the volcanic aerosol cloud, which are summarized and discussed in Sect. 4.

2 Methods and data

2.1 Methods

2.1.1 Experimental protocol

There is a degree of uncertainty over the thickness of the injected SO₂ cloud, based on available measurements. Therefore, different modeling centers may have selected different simulated injection altitudes for the Pinatubo eruption in the past. Within Dhomse et al. (2020), UM-UKCA set the SO₂ injection altitude at 21–23 km based on the altitude of the first detection of the Pinatubo cloud at Mauna Loa (Antuña et al., 2002). Further UM-UKCA analysis by Shallcross (2020) demonstrated improved model correspondence with the July–August 1991 Mauna Loa lidar measurements when running the model with “pre-nudged free-running” rather than the “approximate QBO free-running” (QBO: quasi-biennial oscillation) approach used in Dhomse et al. (2020). Sheng et al. (2015b) performed atmospheric simulations of the Pinatubo eruption with AER 2-D 300 by varying the emission parameters and found agreement with several observations by injecting 14 Tg of SO₂ with a vertical distribution peaking at 18–21 km. Similar emission parameters (10–12 Tg of SO₂ at 18–20 km) were used in Mills et al. (2016) with CESM1-WACCM. Niemeier et al. (2009) showed comparable aerosol optical depth and effective radius with satellite and lidar measurements, simulating with MAECHAM5-HAM the injection of 17 Tg of SO₂ at about 24 km together with 100 Tg of fine ash at about 21 km. Stenchikov et al. (2021) simulated with WRF-Chem v3.7.1 the same amounts of SO₂ and ash but centered at 17 km, showing that the radiative heating of ash can raise the sulfur cloud by 7 km during the first week of the eruption. These differences motivated the design of the ISA-MIP HErSEA intercomparison.

The HErSEA Pinatubo experiment design includes five different emission scenarios considering different amounts and altitudes of injection of SO₂, as summarized in Fig. 1. The first three emission scenarios describe injections at medium altitude (between 21–23 km) of an amount of SO₂ that varies from the lowest value of 5 Tg S (Low-22 km) to a medium value of 7 Tg S (Med-22 km) and the highest value of 10 Tg S (High-22 km). The medium-injection scenario (7 Tg S in the form of SO₂) has three different injection altitude settings: Med-22 km, as discussed; another shallow one at lower altitudes (18–20 km, Med-19 km); and one over a deep altitude range (18–25 km, Med-18–25 km).

The Mt. Pinatubo-like eruption is timed on 15 June 1991. SO₂ is injected in models in a single grid cell close to the Pinatubo location (15^∘ N, 120^∘ E) and at the prescribed altitudes, with the precision given by the specific vertical and horizontal model resolution (Table S1 in the Supplement). UM-UKCA provided an additional set of simulations, called meridional-spread injection simulations, and the EMAC simulation differs from the protocol: this differentiation is highlighted by the addition of a ^* after the model name. In UM-UKCA^*, SO₂ is injected at Mt. Pinatubo longitude and in a latitude range between 0^∘ and 15^∘ N (12 model grid boxes), a common strategy (Dhomse et al., 2014; Mills et al., 2016) to match the initial southward spread of the aerosol cloud (Bluth et al., 1992). In EMAC (we will EMAC^* only in the figures and tables), volcanic SO₂ injections are entered at one single point in time as 3D mixing ratio perturbations derived from satellite data using an inventory for the period 1990 to 2019 (https://doi.org/10.26050/WDCC/SSIRC_3). For the Pinatubo period also the eruptions of Cerro Hudson (10 August 1991), Spurr (27 June 1992), and Lascar (18 April 1993) are included in EMAC. The amount of SO₂ injected is 8.5 and 0.65 Tg S for Pinatubo and Cerro Hudson, respectively, and top heights of the volcanic plumes are approximately 23 and 18 km.

All models are radiatively coupled to the volcanically enhanced stratospheric aerosol in order to resolve the composition–radiation–dynamics interactions. Previous model studies (e.g., Young et al., 1994; Timmreck et al., 1999; Aquila et al., 2012; Sukhodolov et al., 2018) showed that inclusion of the interaction between volcanic sulfate aerosol and radiation is essential for a reliable simulation of the transport of the volcanic cloud. Radiative heating of ash and SO₂ is also important for the initial uplift of the volcanic cloud (Lary et al., 1994; Young et al., 1994; Gerstell et al., 1995), but the contribution of SO₂ is smaller than that of ash, in the first week, or sulfate aerosols, in the subsequent weeks (Stenchikov et al., 2021). About 80 Tg of ash was injected during the Pinatubo eruption (Guo et al., 2004b). However, both ash and SO₂ radiative effects are not included in all model simulations as it is outside the scope of the project, which focuses on the long-term evolution of the Pinatubo volcanic cloud.

Modeling groups performed transient AMIP-type (Atmospheric Model Intercomparison Project) (Gates et al., 1999) runs of the Mt. Pinatubo eruption in which sea surface temperatures and sea ice extent are prescribed as monthly climatologies from the Met Office Hadley Center Observational data set (Rayner et al., 2003). Boundary conditions are also prescribed for greenhouse gases and ozone-depleting substances as recommended for the SPARC CCMI (Stratosphere–troposphere Processes And their Role in Climate Chemistry–Climate Model Initiative) hindcast scenario REFC1SD (Eyring et al., 2013), in order to match those for the time period. The evolution of the quasi-biennial oscillation (QBO) must be consistent through the post-eruption period, as it affects the dispersion of the volcanic plume to mid-latitudes (Trepte and Hitchman, 1992; Baldwin et al., 2001; Punge et al., 2009) and consequently the size distribution and lifetime of stratospheric aerosols (Hommel et al., 2015; Pitari et al., 2016b; Visioni et al., 2017). Accordingly, models with internally generated QBO re-initialized it in order to be consistent with the actual meteorological conditions or used specified dynamics approaches (e.g., Telford et al., 2008). All groups submitted a three-member ensemble for each different injection setting, except for ULAQ-CCM and EMAC, which submitted only one realization. The generation of the ensemble for each model is explained in the respective sections describing the model. Unless otherwise specified, all results shown refer to the ensemble mean.

Figure 1 Graphical representation of injection setting parameters. The reddish boxes represent an injection of 5, 7, and 10 Tg S in the form of SO₂ centered at 22 km; the blue and light-blue boxes represent the injection of 7 Tg S in the form of SO₂ for injection altitudes centered at 19 km and one deep injection between 18 and 25 km.

Cerro Hudson simulations

To evaluate the role of the Cerro Hudson eruption, we performed two additional simulations with the ULAQ-CCM model that, while outside the scope of ISA-MIP, helped clarify some issues raised by the initial results. The two simulations add the Cerro Hudson eruption to the Med-22 km experiment with lower and upper estimates of SO₂ injection based on the Neely III and Schmidt (2016) and MSVOLSO₂L4 inventory (Carn, 2022), respectively. The additional eruption consists of the injection of SO₂ with a uniform vertical distribution on 10 August 1991 in the grid cell corresponding to the Cerro Hudson location (45.9^∘ S, 72.9^∘ W). The lower-end emission, termed Med-22 km + Low-Hud, includes 1.5 Tg of SO₂ between 11 and 16 km, and the upper-end emission Med-22 km + High-Hud includes 4 Tg of SO₂ at 12–18 km.

2.1.2 Participating models

The ISA-MIP multi-model ensemble includes simulations from five global aerosol models: ECHAM6-SALSA, ECHAM5-HAM, SOCOL-AERv2, ULAQ-CCM, UM-UKCA. In addition closely related simulations from a sixth model, EMAC, are considered. The main characteristics of the participating models are reported in Table 1. ECHAM5-HAM, SOCOL-AERv2, and EMAC are based on the same general circulation model (GCM), ECHAM5 (Giorgetta et al., 2006), but with different horizontal and/or vertical resolutions, while ECHAM6-SALSA uses the updated version ECHAM6.3 (Stevens et al., 2013); all have different chemical and aerosol modules.

ECHAM6-SALSA

ECHAM6-SALSA (ECHAM6.3-HAM2.3-MOZ1.0) is an interactive aerosol–chemistry–climate model based on the ECHAM6.3 general circulation model (Stevens et al., 2013). A T63L95 resolution was used in ECHAM6-SALSA simulations, which corresponds to an approximately 1.9${}^{\circ }\times \mathrm{1.9}{}^{\circ }$ horizontal grid and 95 vertical layers reaching up to 80 km. The QBO is internally resolved by the model (Laakso et al., 2022). The GCM is interactively coupled with the HAMMOZ aerosol–chemistry model (Schultz et al., 2018), which is a combination of the Hamburg Aerosol Model (HAM) and the Model for OZone And Related chemical Tracers (MOZART) chemistry model. However MOZART was not used in the simulations for this study, and OH and ozone concentrations were prescribed by a monthly mean climatology; a simplified sulfate chemistry scheme of HAM was used. The aerosol model HAM calculates the emissions, removal, and radiative properties of aerosol. It simulates five major global aerosol compounds: sulfate, organic carbon, black carbon, sea salt, and mineral dust. The aerosol emissions from anthropogenic sources were based on the Community Emission Data System (CEDS) for the CMIP6 anthropogenic emission inventory. Sea salt and dust emissions were calculated online. Aerosol microphysics were calculated by the sectional aerosol module SALSA. A detailed description of the model is given in Kokkola et al. (2018). SALSA describes aerosols using 10 size bins in size space, and the 7 largest bins are separated into externally mixed soluble and insoluble populations. Ensemble members were produced by using insignificantly different values for one of the tuning parameters (the rate of snow formation by aggregation) for January 1991 of each ensemble member.

ECHAM5-HAM

ECHAM5-HAM has the ECHAM5 GCM (Giorgetta et al., 2006), used as a high-top model in the middle atmosphere (MA) version, and is interactively coupled to the aerosol microphysical model HAM (Stier et al., 2005). The horizontal resolution is about 2.8^∘ in longitude and latitude, in a spectral truncation at wave number 42 (T42), with 90 vertical layers up to 0.01 hPa (about 80 km) and an interactive simulation of the QBO. The aerosol microphysical model HAM (Stier et al., 2005) calculates the oxidation of sulfur and sulfate aerosol formation, including nucleation, accumulation, condensation, and coagulation processes. The width of the HAM modes has been adapted to the conditions under a high-sulfur load. The aerosols are prescribed in three modes with a fixed width (Niemeier et al., 2009). HAM was further adopted to stratospheric conditions by applying a simple stratospheric sulfur chemistry above the tropopause (Timmreck, 2001; Hommel et al., 2011). ECHAM prescribes oxidant fields of OH, NO₂, and O₃ on a monthly basis, as well as photolysis rates of OCS, H₂SO₄ SO₂, SO₃, and O₃. The sulfate was radiatively active for both SW and LW radiation and coupled to the radiation scheme of ECHAM. Further details are described in Niemeier et al. (2021). The ensemble members were produced by increasing the stratospheric horizontal diffusion from one level to the next above on 1 January of the year of the eruption. The parameter generating a different dynamical state is perturbed between 1.0, 1.0001, and 1.001.

SOCOL-AERv2

SOCOL-AERv2 is an interactive aerosol–chemistry–climate model that is also based on the ECHAM5 GCM but coupled to the MEZON chemistry (Egorova et al., 2003) and AER sulfate aerosol microphysics (Weisenstein et al., 1997) modules. The model version used here has a horizontal resolution of about 2.8^∘ in longitude and latitude (T42) and 39 vertical layers up to 0.01 hPa. Because of the coarse vertical resolution (∼1.5 km in the lower stratosphere), the QBO is nudged to the observed equatorial wind profiles. The chemistry module calculates the interactions of 89 chemical species of the oxygen, hydrogen, nitrogen, carbon, chlorine, bromine, and sulfur groups in gas-phase, photolysis, and heterogeneous reactions, including reactions in/on aqueous sulfuric acid aerosols. The sulfate aerosol module resolves the aerosol particles in 40 size bins (the highest aerosol size resolution compared to other participating models), ranging in dry radius from 0.39 nm to 3.2 µm, and calculates nucleation, condensational growth, evaporation, coagulation, and sedimentation of sulfate aerosol bins. H₂SO₄ weight percent is calculated online based on actual temperature and relative humidity. Dry and wet deposition of species are interactively calculated based on actual meteorological conditions in the model (Feinberg et al., 2019). Modeled aerosols and chemical species are coupled with the shortwave- and longwave-radiation schemes. Aerosol radiative properties are treated following a lookup-table approach with precalculated values using Mie theory for actual H₂SO₄ weight percent and temperature. All boundary conditions follow the recommendations of ISA-MIP (Timmreck et al., 2018). Three ensemble members were produced by scaling the global CO₂ concentration by ±0.05 %, which started in January 1991 and was maintained for the whole simulation. Besides the 39-level version, SOCOL-AERv2 can also be run on 90 levels, as the other two ECHAM5-based participating models ECHAM5-HAM and EMAC. However, increased resolution more than doubles the computational expenses of the already heavy calculations of interactive chemistry and highly resolved sectional aerosol microphysics. Therefore, the model is mostly used in the 39-level configuration. To test the effects of increased resolution, SOCOL-AERv2 has been additionally used here for the Low-22 km experiment with the 90 levels instead of the 39 reference levels. With this configuration, the model has been spun up to the conditions of 1991. Besides changed resolution, all other settings have been kept the same.

ULAQ-CCM

ULAQ-CCM (University of L'Aquila Chemistry Climate Model) is a global-scale climate–chemistry coupled model with a horizontal resolution of 5${}^{\circ }\times {\mathrm{6}}^{\circ }$ (T21) and 126 log pressure levels (approximate pressure altitude increment of 568 m), from the surface to the mesosphere (0.04 hPa). However, the QBO is not internally resolved and is nudged to observed values (Morgenstern et al., 2017), and its future values are repeated from the historical time series. The chemistry module includes medium- and short-lived species (O_x, NO_y, NO_x, CHO_x, Cl_y, Br_y, SO_x) and the major component of stratospheric and tropospheric aerosols (sulfate, nitrate, organic and black carbon, soil dust, sea salt, polar stratospheric clouds). The microphysical code for aerosol formation and growth includes a gas–particle conversion scheme, homogeneous and heterogeneous nucleation, coagulation, condensation, and evaporation (Pitari et al., 2002, 2016a). It also includes heterogeneous chemical reactions on sulfuric acid aerosols and polar stratospheric cloud particles; both heterogeneous and homogeneous upper-tropospheric formation processes are also included (Visioni et al., 2018a). The aerosol module calculates the aerosol extinction, asymmetry factor, and single-scattering albedo, given the calculated size distribution of the particles for different wavelengths, and they are passed daily to the radiative transfer module, which is a two-stream delta-Eddington approximation model (Toon et al., 1989).

UM-UKCA

UM-UKCA model simulations are performed using the Global Atmosphere 4.0 configuration (Walters et al., 2014, GA4) of the UK Met Office Unified Model (UM v8.4) general circulation model with the UK Chemistry and Aerosol chemistry-aerosol sub-model (UKCA). The GA4 atmosphere model has a horizontal resolution of 1.875${}^{\circ }\times {\mathrm{1.25}}^{\circ }$ and 85 vertical levels (N96L85) ranging from the surface to about 85 km, with an interactive simulation of the QBO. The UM-UKCA configuration adapts GA4 with aerosol radiative effects from the interactive GLOMAP aerosol microphysics scheme and ozone radiative effects from the whole-atmosphere chemistry, which is a combination of the detailed stratospheric chemistry and simplified tropospheric chemistry schemes (Archibald et al., 2020). The GLOMAP stratospheric aerosol microphysics scheme is described in Dhomse et al. (2014), and the model setup is described in Dhomse et al. (2020). Briefly, the model uses the GLOMAP aerosol microphysics module coupled with the troposphere–stratosphere chemistry scheme, and modeled aerosols are coupled with the radiation scheme. The model also uses greenhouse gas (GHG) and ozone-depleting substance (ODS) concentrations from the Ref-C1 scenario used in the CCMI-1 (Morgenstern et al., 2017) activity. Simulations are performed in atmosphere-only mode, and CMIP6-recommended sea surface temperatures and sea ice concentrations that are obtained from https://esgf-node.llnl.gov/projects/cmip6/ (last access: 25 March 2021) are used. Three ensemble members were initialized using the fields of 3 model years of 20-year time-slice simulations prior to 1990 that gave a QBO transition approximately matching that of ERA-Interim reanalysis (for more details, see Dee et al., 2011; Dhomse et al., 2020).

EMAC

EMAC is the ECHAM5 general circulation model coupled with the Modular Earth Submodel System Atmospheric Chemistry (Brühl et al., 2015, 2018). The resolution is T63/L90, i.e., about 1.9^∘ latitude and longitude and 90 layers up to about 80 km with a vertical resolution of about 500 m near the tropopause. The QBO is internally generated but slightly nudged to observations compiled by the Free University of Berlin. Below 100 hPa and above the boundary layer dynamics and temperature are nudged to ERA-Interim. It contains comprehensive gas-phase and heterogeneous chemistry. The applied aerosol module GMXE (Pringle et al., 2010) accounts for seven modes using lognormal size distributions (nucleation, soluble and insoluble Aitken, accumulation, and coarse modes). The boundary between accumulation mode and coarse mode, a model parameter, is set at a dry particle radius of 1.6 µm to avoid too-fast sedimentation of a too-large coarse-mode fraction in case of major volcanic eruptions. Optical properties for the types sulfate, dust, organic carbon and black carbon (OC and BC), sea salt (SS), and aerosol water are calculated using Mie-theory-based lookup tables for each mode consistent with the selected size distribution widths of the modes. This also means that no overall effective radius is used. The resulting total optical depths, single-scattering albedos, and asymmetry factors are used in radiative transfer calculations which feed back to atmospheric dynamics. The results from EMAC were taken from an existing 30-year transient simulation for comparison (Schallock et al., 2021).

Vehkamäki et al. (2002)Vehkamäki et al. (2002)Vehkamäki et al. (2002)Vehkamäki et al. (2002)Pitari et al. (1993)Mann et al. (2010)

Table 1 Main chemical, microphysical, and dynamic characteristics of the participating models.

^* Models with spatially spread SO₂ injections.

2.2 Observation data sets

2.2.1 AVHRR

The Advanced Very High Resolution Radiometer (AVHRR/2) is a space-borne sensor that measures the reflectance of the Earth in five spectral bands covering visible and infrared wavelengths (0.63, 0.86, 3.7, 11, 12 µm). AVHRR/2 instrument was on board the polar-orbiting satellites (POES) NOAA-11 that provided global coverage data with a resolution of 1.1 km and a frequency of Earth scans of twice per day (https://www.avl.class.noaa.gov/release/data_available/avhrr/index.htm, last access: 12 January 2023). The data used here are on a 1${}^{\circ }\times {\mathrm{1}}^{\circ }$ grid as monthly averages (as archived at the NOAA's National Climatic Data Center). As in Long and Stowe (1994) and Aquila et al. (2012), the stratospheric optical depth at 0.5 µm is calculated by removing monthly mean background values (June 1989 to May 1991) from AVHRR observations. The optical depth at 0.5 µm is retrieved through a radiative transfer surface/atmosphere model (RAO et al., 1989); therefore, combined with the previous assumption, AVHRR cannot detect the changes in stratospheric aerosol optical depth (AOD) smaller than 0.01 but can detect values up to 2.0 (Russell et al., 1996).

2.2.2 SAGE II

The Stratospheric Aerosol and Gas Experiment II (SAGE II) is a satellite-based sun photometer that was launched in October 1984 aboard the Earth Radiation Budget Satellite (ERBS) and retired in August 2005. The instrument measures the extinction of the solar radiation through the limb of the Earth’s atmosphere in seven channels ranging from 385 to 1020 nm, with a global coverage from 80^∘ S to 80^∘ N latitude and a vertical resolution of 1 km for the retrieved data (Mauldin et al., 1985). We used the effective radius and the surface area density of aerosol particles from SAGE II version 7.0 (Damadeo et al., 2013; NASA/LARC/SD/ASDC, 2012b). The SAD (and thus the effective radius) is derived by a method that is a linear mix between the Thomason et al. (1997) method, which is valid for the 525–1020 nm extinction ratio below 1.5, and the Thomason and Burton (2008) method for ratios above 2.0 (Damadeo et al., 2013). Both methods assume that aerosols are spherical droplets of H₂SO₄–H₂O solution with a constant composition of 75 % H₂SO₄ and 25 % H₂O by weight. The Thomason et al. (1997) method uses the principal component analysis to derive the SAD from a linear combination of four aerosol extinction measurements (386, 452, 525, 1020 nm). In the Thomason and Burton (2008) method, SAD is derived from the 525 and 1020 nm channels using an empirical parameterization based on the 525–1020 nm extinction ratio.

The stratospheric sulfate burden is taken from the SAGE-3λ data set (ftp://iacftp.ethz.ch/pub_read/luo/CMIP6/, last access: 12 January 2023) that was compiled for Phase 6 of the Coupled Model Intercomparison Project (CMIP6). H₂SO₄ particle number density (and other secondary products not used here) is derived via the SAGE-3λ algorithms that assume a single mode lognormal size distribution of stratospheric aerosol where number density, mode radius, and width are obtained by fitting the SAGE II extinction coefficients at three wavelengths (452, 525, and 1024 nm) (Revell et al., 2017).

2.2.3 HIRS

The High Resolution Infrared Radiation Sounder (HIRS) is an infrared-scanning radiometer that has been onboard several NOAA platforms starting with the first satellite of the Television Infrared Observation Satellite series (TIROS-N), followed by NOAA-6 up to NOAA-19 (Borbas and Menzel, 2021). It measures the reflectance of the Earth in 19 infrared channels (3.7 to 15 µm) and 1 solar channel (0.69 µm) with a spatial resolution at nadir of 20.4 km on HIRS/2. Baran and Foot (1994) used HIRS/2 cloud-cleared radiances at 8.3 µm (NOAA-10/12) and 12.5 µm (NOAA-11) to retrieve the column number density of sulfuric acid aerosols from May 1991 to November 1993. Among the assumption and the approximations, the stratospheric aerosols are assumed to be 75 % H₂SO₄ and 25 % H₂O, with a spectral transmittance based on dustsonde measurements by Deshler et al. (1992) and a single-scattering albedo calculated from Mie theory by integrating the extinction and scattering coefficients over a lognormal size distribution using a mode radius 0.35 µm and a normalized standard deviation of 1.6 (Baran and Foot, 1994). The data cover the latitudes from 80^∘ N to 80^∘ S and all longitudes with 5^∘ of resolution and are affected by a systematic error of 10 % due to the sensitivity of the retrieved method and uncertainties in the background.

2.2.4 OPC

The University of Wyoming balloon-borne Optical Particle Counter (OPC) is a spectrometer that measures the light-extinction cross section of the particles using a broadband incandescent light source, developed by Rosen (1964), providing the particle size and the number concentration. The stratospheric aerosol measurements from 1991 to 2012 are made over Laramie (Wyoming) with the so-called OPC40, which can detect particles throughout the size range 0.1–10.0 µm, distinguished in 8 or 12 channels, depending on the instruments (Deshler, 2003). Here we used the revised data set (UWv2.0; http://www.atmos.uwyo.edu/~deshler/Data/Aer_Meas_Wy_read_me.htm, last access: 12 January 2023) of the OPC measurements (Deshler et al., 2019). Surface area density and volume density are calculated from the size distribution derived from particle size and concentration by fitting the data to a unimodal or bimodal lognormal distribution (depending on the number of measurements and on which of the two minimizes the difference between the calculated and the measured number concentration) (Kovilakam and Deshler, 2015).

2.2.5 GloSSAC

The Global Space-based Stratospheric Aerosol Climatology (GloSSAC) is a global and gap-free data set of zonally averaged optical properties of stratospheric aerosols (focused on aerosol extinction coefficient at 525 and 1020 nm) from 1976–2018. It is mainly based on the Aerosol and Gas Experiment (SAGE) and on the Optical Spectrograph and InfraRed Imager System (OSIRIS) and the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO). Ground, airborne, and balloon-based instruments were used to fill major gaps in the data set (Thomason et al., 2018). Here, we used the updated version v2 (NASA/LARC/SD/ASDC, 2012a) from Kovilakam et al. (2020).

3 Results

The various sets of initial conditions of SO₂ injections result in an aerosol cloud with different optical properties depending on the dispersion of the cloud over time and the size of the aerosols produced.

In the following section, we start by analyzing the AOD and how the models reproduce the measured AOD with different volcanic emission source parameters. Since the amount of attenuation depends on the particle number concentrations and size, we then investigated both the magnitude and distribution of the sulfate burden and the size of the sulfate aerosols.

3.1 Aerosol optical depth

The stratospheric AOD simulated by the different interactive aerosol microphysical models is evaluated by comparing it with satellite observations from AVHRR and GloSSAC (Fig. 2). The AOD is calculated at a wavelength of 550 nm in EMAC, ECHAM5-HAM, ULAQ-CCM, and UM-UKCA; 533 nm in ECHAM6-SALSA; 525 nm in SOCOL-AERv2 and GloSSAC; and 600 nm in AVHRR. Differences between those wavelengths are however negligible. GloSSAC provides zonal values with a latitudinal resolution of 5^∘ and uniform spatio-temporal coverage up to the year 1994. As it is mostly based on SAGE II measurements, the instrument saturates for optical depth of about 0.15; therefore it is less accurate in the center of tropical clouds in the first months after the eruption (Russell et al., 1996). Conversely, AVHRR can only measure stratospheric AOD larger than 0.01. Because of the paucity of data points, “global values” when comparing against AVHRR are calculated between 60^∘ S–60^∘ N.

Figure 2 shows the time evolution of the zonal-mean stratospheric AOD for each model and ensemble mean. It is clear that medium and high injection of SO₂ (Med-22 km and High-22 km, respectively) overestimate the stratospheric AOD in the tropics or/and in the Northern Hemisphere (NH) extratropics compared to both observations. The ability to reproduce the observed values in the Southern Hemisphere (SH) extratropics depends on both the model and the injection parameters. UM-UKCA^* and EMAC, contrary to other models, show more southward transport, probably due to the different injection settings (see Sect. 2.1.1). In UM-UKCA^* the meridional-spread emission (0–15^∘ N) accounts for the initial west-southwestward drift of the volcanic cloud (Bluth et al., 1992), contributing to a more hemispherically symmetric aerosol distribution (Dhomse et al., 2014; Mills et al., 2016; Jones et al., 2017). EMAC used a 3D-plume injection and also included smaller eruptions such as that of Cerro Hudson in the Southern Hemisphere in August 1991 (45.9^∘ S, 72.9^∘ W). The additional injection is a 3D-plume injection of 0.65 Tg S in the form of SO₂, whose maximum in terms of mixing ratio is at 18 km, and differs from the two additional cases performed with ULAQ-CCM (2.1.1.1). In ULAQ-CCM, the Med-22 km+Low-Hud includes a similar amount of SO₂ but at lower altitudes compared to the Cerro Hudson eruption in EMAC, and its effect on the stratospheric burden and AOD is negligible. In contrast, Med-22 km+High-Hud enhances them in the Southern Hemisphere, approaching observation, but only for a few months after the eruption (Fig. S6 in the Supplement).

A quantitative comparison with the observations is shown with the use of Taylor diagrams (see Appendix A) in Fig. 3. Model results are compared for the first year after the eruption with both AVHRR and GloSSAC (first row and second row, respectively) and for the second year only with GloSSAC (third row). Three-member ensembles, when provided, are represented with smaller circles of the same color with respect to the ensemble mean of a specific simulation. In ECHAM6-SALSA, the differences between members of the same scenario are greater than those between scenarios because of differences in local winds at the time of the eruption in each ensemble member. The impact of local winds is weaker when SO₂ is injected over the deep altitude range between 19 and 25 km (blues circles in Fig. 3a and h). There are various sets of initial conditions for SO₂ injections which, depending on the model, are close to the observations. The experiments that best reproduce the observations are those with similar variability to that of the observations, defined by their standard deviations (SDs), higher correlation (COR), and lower root-mean-square difference (RMSD). The values of COR and RMSD for these experiments are summarized in Table 2.

Table 2 Correlation (COR) and root-mean-square difference (RMSD) of the stratospheric AOD calculated between observations and model results for the experiments that best reproduce the observations.

^* highlights models with spatially distributed SO₂ injections.

Figure 2 Time evolution of zonal stratospheric AOD for all models and in Low-22 km (first column), Med-22 km (second column), High-22 km (third column), Med-19 km (fourth column), and Med-18–25 km (fifth column). The last row includes the different scenario simulated by EMAC^* and the two observations used for comparison: GloSSAC and AVHRR. AOD is calculated at a wavelength of 550 nm in ECHAM5-HAM, EMAC, ULAQ-CCM, and UM-UKCA; 533 nm in ECHAM6-SALSA; 525 nm in SOCOL-AERv2; 525 nm in GloSSAC; and 600 nm in AVHRR. ^* Models with spatially spread SO₂ injections.

Figure 3 Taylor diagrams for the global stratospheric AOD. Zonal monthly mean values for different time periods have been used to calculate the standard deviation, correlation, and centered root-mean-square difference between model experiments and measurements. In the first row, model results are compared with respect to AVHRR over the period June 1991 to May 1992, in the second row with respect to GloSSAC over the period June 1991 to May 1992, and in the third row with respect to GloSSAC over the period June 1992 to May 1993 (See Appendix A1 for more details). ^* Models with spatially spread SO₂ injections.

During the first year after the eruption, all models show better agreement with AVHRR than GloSSAC: correlations range between 0.73 and 0.78 with AVHRR versus 0.54 and 0.82 with GloSSAC, for which RMSDs are also higher. In ECHAM6-SALSA, SOCOL-AERv2, and ULAQ-CCM, the injection of 7 Tg S in the form of SO₂ closer to the tropopause is a good compromise between the too-high and too-low stratospheric AOD produced in the tropics by an injection of 5 and 10 Tg S in the form of SO₂, respectively, and this scenario also produces a better southward and northward transport (Fig. 2). The best set of initial parameters also depends on the observation considered for comparison: in ECHAM6-SALSA Med-18–25 km and Med-19 km reproduce better AVHRR and GloSSAC measurements, respectively, and in the comparison with GloSSAC the correlation increases, and the RMSD decreases over time (Fig. 2a5). For SOCOL-AERv2 and ULAQ-CCM, Med-19 km is in good agreement with both AVHRR and GloSSAC in the two different periods considered (Fig. 2c4 and d4). During the first year after the eruption, the correlation between Med-19 km and the observations is higher for ULAQ-CCM (0.84 and 0.74 compared with AVHRR and GloSSAC, respectively) as it better reproduces the tropical confinement, while in the following year (June 1992–July 1993), in SOCOL-AERv2 comparable values of stratospheric AOD persist for longer in the extratropics compared with GloSSAC (correlation of 0.86). In ECHAM5-HAM the injection at 21–23 km results in a comparable stratospheric AOD in the tropics and SH extratropics compared to both observations but overestimates northern-hemispheric (NH) extratropics values by up to a factor of 2 (Fig. 2b1, b2, and b3). The amount of SO₂ to obtain the highest correlation between modeling experiments and observations depends on the observation and on the period considered: High-22 km and Low-22 km when compared with AVHRR and GloSSAC during the first year after the eruption, respectively, Med-22 km when compared with GloSSAC the following year. In UM-UKCA, the point injection and meridional-spread emission agree that Low-22 km better reproduces the stratospheric AOD of both observations during the first year after the eruption, as it shows a good tropical confinement and comparable values in the NH, and for the meridional-spread emission also in the SH (Fig. 2e1 and f1). Therefore, the correlation is higher, and the RMSD is lower for the meridional-spread emission experiment. The poleward transport, especially in the NH, is enhanced in Med-19 km (Fig. 2e4 and f4) and found to have a higher correlation with GloSSAC 1 year after the eruption (COR of 0.86 and 0.47 for UM-UKCA^* and UM-UKCA, respectively). During the first year after the eruption, EMAC has comparable values in the tropics and northern mid-latitudes with respect to AVHRR, while in the southern mid-latitudes the stratospheric AOD is up to twice as large and results in a correlation of 0.79. The correlation decreases to 0.63 when comparing with GloSSAC during the following year because of the more rapid decline in the stratospheric volcanic cloud.

Figure 4 Time evolution of monthly values of the normalized global stratospheric AOD for models (colored lines) and AVHRR and GloSSAC observations (black lines). The dashed gray line represents the $\mathrm{1}/e$ value. The experiments shown are Med-19 km for ECHAM6-SALSA, SOCOL-AERv2, ULAQ-CCM, UM-UKCA, and UM-UKCA^* and Med-22 km for ECHAM5-HAM. For EMAC^*, it refers to the only experiment provided. ^* Models with spatially spread SO₂ injections.

The persistence of the volcanic aerosol in the stratosphere is shown in Fig. 4, which represents the global normalized stratospheric optical depth, calculated as explained at the beginning of Sect. 3.1. The Med-19 km experiment is shown for all models, as it is the experiment which best reproduces the GloSSAC observations after June 1992 for all models, with the exception of Med-22 km for ECHAM5-HAM and EMAC with the only experiment provided. The e-folding time, calculated as the time between the maximum and the 1/e value, is 13 months in AVHRR and 15 months in GloSSAC. This range includes ULAQ-CCM and UM-UKCA, with an e-folding time of 14 months, and UM-UKCA^*, with an e-folding time of 15 months. Lower values were found for SOCOL-AERv2 with 12 months, ECHAM6-SALSA and ECHAM5-HAM with 11 months, and EMAC with 10 months.

3.2 Sulfate burden

Figure 5 shows the time evolution of the global and tropical stratospheric sulfate burden of different injection setups for each model. The results of each model are compared with satellite measurements from HIRS and the SAGE-3λ data set. Large differences are evident in the temporal evolution of the sulfate burden between the aerosol model simulation on one hand and the satellite data set on the other, which show similar values and a similar temporal evolution for the sulfate burden.

Figure 5 Time evolution of monthly values of global and tropical stratospheric sulfate burden in teragrams of sulfur (first and second column, respectively). Each panel refers to the respective model in which the different results of the experiments (colored lines; different line styles for different experiments; see legend on the left) are compared with the HIRS and SAGE-3λ data sets (black lines; see legend on the right). ^* Models with spatially spread SO₂ injections.

In the 6 months following the eruption (July–December, termed the build-up phase), ECHAM6-SALSA, ECHAM5-HAM, SOCOL-AERv2, and ULAQ-CCM best match the global stratospheric sulfate burden of HIRS and SAGE-3λ with the injection 5 Tg S in the form of SO₂ (Low-22 km), a lower amount compared to the one required for a comparable stratospheric aerosol optical depth (Fig. 5a, b, d, and e). For SOCOL-AERv2, Med-19 km also shows values within the uncertainties in the HIRS measurements. However, Low-22 km, and also Med-19 km for SOCOL-AERv2, anticipates the peak and underestimates the tropical burden in ECHAM6-SALSA, ECHAM5-HAM, and SOCOL-AERv2, while the peak is reached later, and larger values are produced in ULAQ-CCM (Fig. 5h, i, k, and l). In UM-UKCA, point and meridional-spread injection show similar results for the global stratospheric sulfate burden and agree with observations with Med-19 km and Med-18–25 km experiments (Fig. 5f and g). The differences between the two strategies emerge in the tropics, where values are lower for point injection experiments due to the lack of aerosols transported to the southern tropics and that are therefore confined to the Northern Hemisphere. For the point injection, Low-22 km and Med-18–25 km approaches SAGE-3λ for the first months and HIRS for the last 3 months of the build-up phase. All the experiments with larger amounts of injected SO₂, including the EMAC experiment with 8.5 Tg S in the form of SO₂, overestimate the measured global sulfate burden; all experiments in ULAQ-CCM and the single scenario in EMAC overestimate the tropical burden, while in ECHAM6-SALSA, ECHAM5-HAM, and SOCOL-AERv2 they overestimate the burden in the NH extratropics (Fig. S5).

In the build-up phase, SAGE-3λ assumes the lowest values and slowly reaches a peak of 5.0 Tg S in December, compared to 5.4 Tg S of HIRS in September. Lower values in SAGE-3λ are related to the saturation effects of the limb-occultation instrument; therefore HIRS measurements are to be considered more reliable for this initial period (Sukhodolov et al., 2018). For EMAC, the injection of 8.5 Tg S in the form of SO₂ produces a sulfate aerosol cloud that peaks in September at 7.0 Tg S, a value comparable to the results of the Med-22 km experiment (performed by the other models), in which 7 Tg S in the form of SO₂ is injected. For SOCOL-AERv2 and UM-UKCA with both injection strategies, Med-19 km shows the best agreement with HIRS in terms of peak and timing of the peak (September for SOCOL-AER, October for UM-UKCA), whereas in Low-22 km and the other experiments it is reached 1 month later. This is followed by ECHAM6-SALSA in October (November only in High-22 km) and ULAQ-CCM in November. ECHAM5-HAM is more sensitive to the altitude of injection: it peaks between October in Med-19 km, November in Med-18–25 km, and December in the experiments with the same altitude of injection (Low-21 km, Med-21 km, and High-21 km); the values of the peak are 14.3 % lower in Med-19 km and 7.1 % lower in Med-18–25 km compared to Med-22 km.

The sensitivity to injection altitude depends on the model: during the build-up phase, the Med-18–25 km and Med-22 km curves coincide in ECHAM6-SALSA and SOCOL-AERv2, and, compared to these experiments, the values in Med-19 km are up to 9 % and 20 % smaller for each model, respectively. In ULAQ-CCM, ECHAM5-HAM, and UM-UKCA, the more SO₂ is injected at lower altitudes the smaller the value of the peak is, but for ULAQ-CCM the peak is only 1 % and 6 % lower in Med-18–25 km and Med-19 km compared to Med-22 km. The value and time of the peak for all models and experiments are summarized in Table S2. In general, when the amount of SO₂ injected is exclusively in the lowest levels or in some vertical levels that include the lowest levels (Med-19 km and Med-18–25 km, respectively), the sulfate burden is lower, and therefore this effect is less pronounced at Med-18–25 km, as the aerosol distribution is more dependent on the balance between gravitational sedimentation in the lower stratosphere and the strength of vertical transport by the Brewer–Dobson circulation, as well as the height of the tropopause.

Figure 6 Time evolution of global stratospheric sulfate burden normalized to the amount of injected SO₂. Each panel refers to the respective model in which the different experiments are compared.

Differences among models and experiments in terms of amount and timing during the build-up phase are influenced by the oxidation of SO₂ by OH that determines the timescale for aerosol formation (Clyne et al., 2021). For this reason, we distinguish between models with prescribed OH (ECHAM6-SALSA and ECHAM5-HAM) and those with interactive OH (SOCOL-AERv2, ULAQ-CCM, UM-UKCA) when looking at the SO₂ evolution. The global normalized SO₂ burden curves (Fig. S4a) coincide for all models with prescribed OH. An exception is Med-19 km in ECHAM6-SALSA, which has lower values and might depend on an early removal through tropopause flux, facilitated by injection near the tropopause. In ULAQ-CCM and UM-UKCA, when comparing High-22 km with Low-22 km we find that a higher injected SO₂ mass produces a longer initial e-folding time for SO₂. The same applies when comparing injections concentrated in a few kilometers (Med-22 km and Med-19 km), i.e., where SO₂ oxidation depletes OH more quickly (Mills et al., 2017), with those where the same amount of SO₂ is injected over a wider altitude band. Consequently, initial values of the stratospheric sulfate burden in Med-18–25 km are slightly higher compared to Med-22 km and Med-19 km.

In order to better understand the models' sensitivity to the different emission scenarios and eventual non-linearities, in Fig. 6 we normalize the resulting global sulfate burden by the amount of SO₂ injected. Thus, in the build-up phase we would expect all the curves for all experiments to reach a value of 1, since no SO₂ and sulfate aerosols have yet been removed from the atmosphere. This will highlight the differences in the aerosol removal (wet removal, deposition, sedimentation) depending on the injection altitude and differences in microphysical growth, especially in the descending phase. Not all models and experiments, however, reach the value of 1: ECHAM5-HAM in Med-19 km and Med-18–25 km and ULAQ-CCM in Med-19 km do not, nor do any experiments in ECHAM6-SALSA, SOCOL-AERv2, and UM-UKCA. This is due to the use of monthly averages for our analyses and the faster removal, near the tropopause, of sulfate aerosol and SO₂ not yet converted to aerosols, especially in Med-19 km and Med-18–25 km experiments. To confirm this, we observe that this is particularly evident in Med-19 km with the lowest injection height. The curves of the experiments with injection between 21–23 km coincide in the build-up phase and the differences emerge later, after 1992: the aerosol lifetime decreases with increasing mass of SO₂ injected (Table S2), which corresponds to the increase in the aerosol size in all models. In UM-UKCA, the lifetime is increased by 1 to 2 months for the meridional-spread emission compared with the point injection. In ECHAM6-SALSA the lifetime increases when increasing the injected SO₂ mass. However, Figs. 3 and S1 show that the differences in results between ensemble members of the same scenarios are larger in ECHAM6-SALSA than in other models. This indicates that differences in aerosol lifetimes between Low-22 km, Med-22 km and High-2 km scenarios are probably not statistically significant in ECHAM6-SALSA. Figure S11a shows the sulfate burden from SOCOL-AERv2 for the Low-22 km experiment calculated with two vertical model resolutions. This figure further confirms the faster removal of volcanic sulfur during the first months after the eruption in SOCOL-AERv2 even in the 22 km injection experiments. The lower-vertical-resolution version shows a much lower burden peak already in late 1991, while the higher-resolution version peaks at exactly the emitted amount of 5 Tg S plus the background value of ∼0.17 Tg S and maintains this peak till early 1992. This is an effect of increased vertical diffusion in the lower-resolution version, which quickly redistributes the volcanic cloud vertically in both directions. This brings some of the volcanic sulfur mass closer to the tropopause and the shallow branch of the Brewer–Dobson circulation, reducing its confinement in the tropical reservoir and enhancing removal from the stratosphere (Brodowsky et al., 2021). This agrees with the results of 22 km experiments of high-resolution ECHAM5-HAM, which also maintain the emitted amount for some months after the eruption (Fig. 6).

Among all models and experiments, the shortest e-folding time of the global stratospheric sulfate burden is 8 months for EMAC; ranges between 10 and 14 months for ECHAM6-SALSA, ECHAM5-HAM, SOCOL-AERv2, and ULAQ-CCM; and reaches the highest values for UM-UKCA, with values between 17 and 23 months, which more closely matches those of HIRS and SAGE-3λ of 21 and 20 months, respectively. The e-folding time of the tropical stratospheric sulfate burden is 12 and 13 months in HIRS and SAGE-3λ and half for the models, with the exception of ECHAM5-HAM for Low-22 km, Med-22 km, and Med-18–25 km, with a longer duration of 9 months, and UM-UKCA, for which it varies between 8 and 14 months, based on the experiments and injection strategy. No model except UKCA can reproduce the observed slow-descent phase during 1992 of the stratospheric sulfate burden, and only the High-22 km scenario approaches the measured values at the end of 1992 for these models, while strongly overshooting them in the preceding months.

Overall, we find that Low-22 km and High-22 km are the experiments that, in all models, better reproduce the observations in the build-up and descent phase, respectively (Figs. 5, S6). The spatio-temporal development of the sulfate burden (Fig. S6) reflects in general that of the AOD (Figs. 2, 3). In the SH, the stratospheric burden shown in SAGE-3λ is not reproduced by the models in Low-22 km; therefore more SO₂ (High-22 km) must be injected for the aerosol cloud to persist for as long as in SAGE-3λ and reach the same values. This way, however, the burden in the NH is overestimated (Fig. S5). There are clear differences in the position of the stratospheric AOD peak, which lies between 5–20^∘ N in the models but around 5^∘ S–10^∘ N in the observations, pointing to differences in the meridional transport in the early phase after the eruption (Fig. 2). In addition, Fig. S11b–c illustrate that the volcanic aerosol mass redistribution between the hemispheres could also be affected by the vertical resolution of the models because it affects the timings of tropical confinement and across-tropopause removal.

Figure 7 Time evolution of the latitudinal partition of the stratospheric sulfate burden. The aerosol mass fraction is calculated with respect to the total burden, for the tropical burden (20^∘ N–20^∘ S) (a, d g), for the burden integrated over the northern mid-latitudes (35^∘–60^∘ N) (b, e, h), and for the burden integrated over the southern mid-latitudes (35^∘–60^∘ S) (c, f, i). The first row includes the experiments with different amounts of SO₂ injected, the second row experiments with different injection altitudes. The third shows the percentage change in the latitudinal partition for all model experiments compared to SAGE-3λ. Experiments are identified here with different line styles; the different colors refer to the models. ^* Models with spatially spread SO₂ injections.

In order to discuss the meridional transport, Fig. 7 shows the aerosol mass fraction of the simulated sulfate burden in the tropics (20^∘ N–20^∘ S), in the northern mid-latitudes (35^∘–60^∘ N), and in the southern mid-latitudes (35^∘–60^∘ S) with respect to the global value, for SAGE-3λ (black line), and for all models and scenarios (first row for the different injection amounts, second row for the different injection altitudes). Tropical confinement (Fig. 7a and d), as shown in the observations, is not captured by ECHAM6-SALSA, ECHAM5-HAM, SOCOL-AERv2, and EMAC, which underestimate the tropical aerosol mass fraction, resulting in a stronger transport to the NH for the first three models and to the SH for EMAC. ULAQ-CCM overestimates the fraction during the first 6 months after the eruption and becomes comparable thereafter. UM-UKCA shows tropical confinement comparable to that of SAGE-3λ for the 21–23 km injection experiments for point injection and shallow and deep injection for meridional-spread emission, otherwise underestimated or overestimated in the other experiments, respectively. However, the similarity between observations and the 21–23 km injection experiments for the UM-UKCA point injection masks the lack of aerosols in the southern tropics (0–20^∘ S) and a higher load in the northern extratropics (0–20^∘ N). Indeed, the fraction of burden for the NH mid-latitudes (Fig. 7b and e) is overestimated, with differences of up to 20 % compared to SAGE-3λ (Fig. 7h), while for the SH (Fig. 7c and f) it is underestimated but to a smaller extent, with differences of 10 % compared to SAGE-3λ (Fig. 7i). The same happens for ECHAM6-SALSA, ECHAM5-HAM, and SOCOL-AERv2. Overall, NH transport is favored in all models at the expense of tropical confinement.

In most models, varying the injected SO₂ mass does not affect the fraction of aerosols transported out of the tropics towards both hemispheres (Fig. 7a, b, and c). The only exception is ECHAM6-SALSA, where an increased injected SO₂ mass increases the tropical confinement, especially in the first 6 months after the eruption. All models, except ULAQ-CCM, show that the tropical confinement is reduced in favor of transport towards both hemispheres when SO₂ is injected below 20 km (Med-19 km). Compared to high-altitude injection settings (>20 km), Med-19 km has the greatest transport in SH. The increase in altitude of injection (Med-22 km and Med-18–25 km) produces a higher confinement in the tropics with a consequent reduced transport toward both hemispheres in ECHAM6-SALSA, SOCOL-AERv2, and UM-UKCA. In ECHAM5-HAM, the strongest confinement is achieved in Med-22 km, while Med-18–25 km shows a similar behavior to Med-19 km as most of the sulfate aerosols are found below 20 km. In ULAQ-CCM differences among the injection settings emerge 6 months after the eruption, and the injection at lower altitudes (Med-19 km) shows a more efficient polewards transport, especially towards the NH.

Figure 8 Time evolution of stratospheric effective radius (µm) in the tropics (a–g) and over Laramie (41^∘ N, 105^∘ W) (h–n). In the panels of the first row, the stratospheric effective radius of the models is calculated between 21–27 km (50–20 hPa) to be compared with the available SAGE II observations. In the panels of the second row, it is calculated between 14–30 km (130–10 hPa) to be compared with the OPC observations. ^* Models with spatially spread SO₂ injections.

3.3 Effective radius and surface area density

Figure 8 shows the time evolution of the observed and simulated stratospheric effective radius in the tropics (20^∘ S–20^∘ N) and over Laramie (41^∘ N–105^∘ W) (calculation of the effective radius and error bar in Appendix A2). In the tropics (Fig. 8a–g) the stratospheric effective radius is calculated as the SAD-weighted average between 21–27 km because of a paucity of tropical measurements below 21 km in SAGE II. Over Laramie (Fig. 8h–n), the stratospheric effective radius is defined as the SAD-weighted average between 14–30 km in order to compare it with in situ OPC measurements (Deshler et al., 2019). Model results are calculated as the value of the nearest grid cell to Laramie; therefore, the ability to reproduce the OPC measurements is more influenced by atmospheric circulation patterns as zonal-mean comparisons discussed earlier and depends also on the horizontal resolution (see Table S1).

Before the eruption, the simulated evolution of the tropical-mean effective radius in most models is almost steady compared to SAGE II. Only ULAQ-CCM reproduces the observed seasonal variation and matches the pre-eruption measurements, resulting in particles with a radius of 0.27 µm, similar to SAGE II (calculated over the 5 months before the eruption). The other models have smaller background particles with a constant value of 0.14 in ECHAM6-SALSA, 0.17 in ECHAM5-HAM, 0.17 in EMAC, 0.15 in SOCOL-AERv2, and 0.10 in UM-UKCA. Over Laramie, ECHAM6-SALSA, ECHAM5-HAM, EMAC, and SOCOL-AERv2 have comparable radii to the OPC ones, while ULAQ-CCM and UM-UKCA lie outside the uncertainty range with larger and smaller radii, respectively. The causes of these differences are unclear; however, an in-depth exploration of the background behavior is out of the scope of this paper and needs to be addressed by studies specifically designed to study aerosol microphysics and transport under volcanically quiescent conditions such as the ISA-MIP background experiment (Timmreck et al., 2018).

After the eruption, all models are able to capture the same decay rate as the SAGE II measurements, remaining flat around the peak reached approximately after October 1991. Most produce a comparable tropical effective radius for about a couple of years, based on different injection settings. The models agree that particle size increases with increasing injected SO₂ mass, with differences from the medium-injection scenario within 15 % in ECHAM6-SALSA and 10 % in ECHAM5-HAM, SOCOL-AERv2, ULAQ-CCM and UM-UKCA. The differences are larger when comparing different injection altitude scenarios, and the corresponding increase in the particle size is model-dependent. In ECHAM6-SALSA and SOCOL-AERv2, High-22 km shows a tropical stratospheric effective radius within 10 % of SAGE II until the end of 1993, peaking, respectively, at 0.47 and 0.49 µm compared to 0.51 in SAGE II. In ECHAM5-HAM, all experiments except High-22 km, which best fits the observed AOD (see Sect. 3.1), produce similar effective radii, ranging between 0.46 and 0.51 µm, and are comparable with SAGE II until the end of 1992. High-22 km differs by larger radii reaching a maximum of 0.56 µm. One year after the eruption, the differences among the different ECHAM5-HAM experiments disappear, and the effective radius decreases more rapidly than in SAGE II. EMAC peaks at 0.33 µm in October, and radii stay around 0.30 µm for less than 1 year. The low bias hides the faster decrease in the effective radius at about 22 km altitude than in most other models, while in the stratosphere below it is similar to observations. In ULAQ-CCM, the effective radius of Med-19 km reproduces the SAGE II measurements with a similar time decrease, as differences stay within 10 % until the end of 1995, while other experiments produce larger particles, with peaks ranging between 0.53 and 0.71 µm. In UM-UKCA, the growth of the effective radius is slower compared to other models, particularly for point injection, but both injection strategies show the slowest decay, which is closest to that of SAGE II. After peaking at different times, the radii between the two injection strategies are similar and range between the smallest value of 0.10 for Med-19 km and the largest value of 0.49 in High-22 km, which is comparable with the observations.

Over Laramie, all experiments of ECHAM6-SALSA, SOCOL-AERv2, and UM-UKCA produce radii within the estimated uncertainties in the OPC measurements for all 5 years in the first two models and after the end of 1991 in UM-UKCA. ECHAM5-HAM and EMAC show comparable values during the pre-eruption phase, but in ECHAM5-HAM radii rise faster compared to the observation during the build-up phase, while in EMAC, after reaching a peak that is about 30 % smaller than that of OPC, the radii assume the smallest values, below the uncertainty. In ULAQ-CCM, all experiments overestimate OPC measurements until early 1992; in particular Med-19 km peaks at 0.78 µm in November 1991, and the effective radius remains at the upper extreme of measurement uncertainty from there on. Increased vertical resolution calculations with SOCOL-AERv2 reveal no difference to the aerosol size before and 1.5 years after the eruption compared to the reference configuration (Fig. 11f–g). During the period of the tropical residence, however, the effective radius noticeably increases due to more aerosol staying in the tropics and the stratosphere and thus available for coagulational growth.

Figure 9 Vertical profile of the effective radius in micrometers (a, d), surface area density (SAD) in square micrometers per cubic centimeter (b, e), and extinction at 0.5 µm in km⁻¹ (c) in the tropics (a–c) and over Laramie (d–e) for Med-22 km in December 1991. Model results are compared with SAGE II and GloSSAC in the tropics and with OPC over Laramie. ^* Models with spatially spread SO₂ injections.

Figure 9 summarizes the information regarding the vertical distribution of the effective radius, SAD, and extinction at 0.5 µm for the Med-22 km experiment, in the tropical area (20^∘ S–20^∘ N), and over Laramie, 6 months after the eruption. A corresponding figure including all available experiments is shown in Fig. S10. By looking at the vertical profiles of various quantities, biases that are hidden in integrated variables emerge. Figure 9c reveals that the vertical profiles differ not only between models and observations but also strongly between the observations themselves.

In the tropics, the effective radius peaks between 100–50 hPa in ECHAM6-SALSA, EMAC, and ULAQ-CCM and between 50–20 hPa in ECHAM5-HAM and UM-UKCA as in SAGE II, with values within 30 % of that measured, except for ULAQ-CCM, where the radii are up to 4 times larger. In UM-UKCA, the peak of SAD for point injection is centered at higher altitude, around 30 hPa compared to 20 hPa for meridional-spread emission, and with smaller values. SOCOL-AERv2 shows good agreement with SAGE II between 100–20 hPa, with values that remain constant around 0.44 µm above 70 hPa. The tropical SAD simulated by the models follows the same vertical distribution as that of SAGE II, and all models have a peak between 50–20 hPa, with the exception of EMAC, whose peak is around 50 hPa. In that range of altitudes, the values of the SAD are comparable with the observations for SOCOL-AERv2 and ULAQ-CCM for most of the attitudes and are up to 2 times larger in the other models.

The tropical extinction follows the same distribution of the SAD. In this case, the extinction is compared with SAGE II and GloSSAC, and large differences exist between them: below 20 hPa the extinction in GloSSAC is larger than in SAGE II, and the differences increase with decreasing height up to 100 % compared to SAGE II because of its gap-filling with ground-based measurements (Thomason et al., 2018; Kovilakam et al., 2020). Above 70 hPa, around the lower bound of the injection altitude, models' extinction is even larger than GloSSAC: ECHAM6-SALSA, SOCOL-AERv2, and ULAQ-CCM approach the measurements at the limit of maximum uncertainty around at 70–25 hPa, and EMAC does so between 40–20 hPa, while ECHAM5-HAM and UM-UKCA overestimate measurements up to twice their value. Below 70 hPa, all models underestimate the GloSSAC data, but the models' extinction is still larger than that of SAGE II, with the exception of EMAC, which shows the greatest extinction below 50 hPa, where it peaks. Considering that the SAD depends on the size and the number of particles, we can assume, for the models that show a comparable radius and a larger SAD compared to SAGE II in the tropics, that they overestimate the number of optically active particles and therefore show a larger extinction (ECHAM5-HAM and UM-UKCA).

Over Laramie, the vertical distribution of the effective radius is within the error bar of the OPC measurements up to 20 hPa in ECHAM6-SALSA, ECHAM5-HAM, and SOCOL-AERv2, while ULAQ-CCM produces larger particles, especially below 50 hPa. In EMAC the effective radius is at the lower limit of the uncertainty but is the only model able to reproduce the vertical profile of the SAD from OPC measurements for most of the altitudes. The models that showed faster transport in the northern mid-latitudes overestimate the observed SAD for most of the altitudes.

The ability to reproduce the observations also depends on the period considered (Figs. S8 and S9): in the first months after the eruption, models and observations show large differences, especially for SAD and extinction, which are overestimated at both latitudes considered. This may be related both to the sensitivity to the actual meteorological conditions that climate models are unable to accurately replicate and to the absence in HErSEA simulations of volcanic ash injection that could remove some of the initial SO₂ gas or affect the local winds and the SO₂ dispersion Zhu et al. (2020); Dhomse et al. (2020); Kloss et al. (2021); Niemeier et al. (2021). This sensitivity to the initial conditions of SO₂ injections decreases the more time passes after the eruption. One year after the eruption, the models still show a vertical profile of the effective radius comparable to observations, while the simulated SAD starts to decrease everywhere after 6 months from the eruption, underestimating tropical values but still overestimating OPC measurements.

4 Discussion

With the use of Taylor diagrams, we highlighted the experiments that better match the observations in terms of stratospheric AOD, in two different time periods, based on the reliability of the measurements. Each model requires different injection scenarios to reproduce the observations, due to differences in the transport and microphysical processes and their mutual interaction. Even considering the best set of initial parameters based on AOD (Fig. 2), differences with observations more or less persist in the models, and we cannot unequivocally define a “best” model as that varies depending on the variable considered and the timing of the observation.

Comparing the results of the models between the experiments with the same injection setup, we observe a large difference between models in reproducing the stratospheric optical depth compared to the similar evolution of the global stratospheric sulfate burden. It is hard to disentangle the transport and the microphysics contribution to the differences in the considered variables, i.e., what fraction of it depends on microphysical schemes or different dispersion of the aerosol cloud. We first considered the contribution of SO₂ oxidation by OH to differences in the timing of the peak for the stratospheric sulfate burden (Fig. 5) and, consequently, AOD (Fig. S2). For models with prescribed OH, differences in the stratospheric rate of SO₂ conversion may depend on the injection altitude, due to an earlier removal through the tropopause flux when the injection is closer to the tropopause. For models with interactive OH we observe a longer e-folding time for higher mass of SO₂ injected and when injected in a narrow altitude range (Med-22 km vs. Med-18–25 km). Due to the availability of only monthly values, some observations of the SO₂ behavior at a more finely resolved temporal scale are not possible here. Furthermore, since the lifetime of sulfate depends on OH concentration and transport and mixing into adjacent grid boxes, when comparing different models, the timing of the peak cannot be simply related to the treatment of OH.

However, we find a common problem in transport, either too fast from the tropics to high northern latitudes (ECHAM6-SALSA, ECHAM5-HAM, SOCOL-AERv2), confined in the NH (UM-UKCA for point injection), or too confined to the tropics (ULAQ-CCM). The different tropical confinement can be affected by a different vertical advection scheme between ULAQ-CCM and the other models, based on the same dynamical core ECHAM5 or ECHAM6. Here, the tropical confinement depends on the different horizontal resolution (Niemeier et al., 2020), while the particular definition of the tropical pipe (see Waugh et al., 2018) may also strongly affect this conclusion. The vertical resolution of a model can also affect the transport from the tropics to high northern latitudes: Brodowsky et al. (2021) showed for the SOCOL-AER model that a longer tropical confinement was found with increased vertical resolution. Hence, the transport to NH and SH can depend on model version and injection setting: the previous MAECHAM5-HAM simulation of the Mt. Pinatubo eruption by Niemeier et al. (2009) showed a similar pattern for the stratospheric AOD compared to AVHRR and SAGE II, by injecting a mid-range amount of SO₂ between the Med-22 km and High-22 km experiments (8.5 Tg S) into one grid box at the location of Pinatubo and a model layer around 24 km, but assuming fewer vertical levels without internally generated QBO. The Typhoon Yunya, which cannot be reproduced with coarse resolution in models, might have played a role in the equatorward transport of the volcanic cloud as well, causing a stronger transport into the SH than in most model results. Better transport to the SH showed EMAC, which has been nudged to the real meteorological conditions and the UM-UKCA version with emissions between 15^∘ N and the Equator.

The meridional transport in the models depends on the vertical wind structure and on the vertical distribution of the simulated volcanic cloud in the first months after the eruption. Labitzke and McCormick (1992), based on SAGE II measurements, showed for the early post-Pinatubo period an upper transport regime (above 20 km) in which aerosols remain confined to the tropical reservoir spreading between 30^∘ N and 10^∘ S and a lower transport regime (below 20 km) in which aerosols mainly spread to northern high latitudes. Between August and September, aerosols above 20 km spread across most of the SH, reaching latitudes of 50^∘ S, followed in November and December by an enhancement in the NH due to the transition from boreal summer to winter circulation in the middle and upper stratosphere. Most of the models show that a faster transport in the NH is favored when aerosols are mainly distributed in the lower transport regime (Timmreck et al., 1999). The lower-stratospheric part of the injection profiles is also strongly affected by the inconsistencies between the modeled and real tropopause heights at the time of eruption (Brodowsky et al., 2021). This effect can be additionally enhanced in the models with low vertical resolution (Fig. S11). We note that the strength of the meridional transport is also seasonally dependent, and therefore eruptions happening in other seasons would result in different distributions of the aerosol cloud (Visioni et al., 2019; Toohey et al., 2011). We find that the injection rate does not affect the fraction of aerosols transported out of the tropics towards both hemispheres, with the exception of ECHAM6-SALSA, where an increased injected SO₂ mass increases the tropical confinement, especially in the first 6 months after the eruption. This is probably due to a stronger radiative interaction from the absorption of more longwave radiation by larger particles. The behavior of the other models is consistent with the findings of Young et al. (1994) and Aquila et al. (2012), where the aerosol heating by absorption of the infrared radiation induces a lofting and a divergent motion that affects only the initial transport (within 1 month) of the aerosols towards and within both the northern and southern tropics.

Even when models and measurements look comparable for the integrated variables (Figs. 8 and S2), these similarities hide the models' inability to reproduce the observed vertical structure depending on the latitude and time period after the eruption under consideration (Figs. 9, S8, and S9). Most models take up to 6 months before they can reproduce the vertical structure of effective radius, SAD, and extinction in the tropics and up to a year at mid-latitudes. The vertical distribution of SAD and effective radius in three moments identifying the build-up, maximum, and descent phase of the evolution of the sulfate burden (September and December 1991 and June 1992, respectively) shows an initial overestimation of the observations and an underestimation 1 year after the eruption. The lack of ash co-emission, a process not included in HErSEA simulations, could be crucial in the first days/month to better reproduce the initial cloud evolution (Stenchikov et al., 2021). On one hand, the ash may have removed parts of the initial sulfur cloud through the SO₂ or H₂SO₄ uptake on these coarse particles, which have a significant fall velocity (Zhu et al., 2020); on the other hand, the presence of smaller ash particles causes greater heating and vertical lofting of the volcanic cloud (Niemeier et al., 2021; Kloss et al., 2021), which could result in slower meridional transport and longer lifetimes of stratospheric volcanic aerosols, depending on the latitude and injection altitude of SO₂ (Niemeier et al., 2009; Stenchikov et al., 2021). Aberystwyth lidar measurements from Vaughan et al. (1994) show a signature of depolarizing particles around 16 km between November and December 1991. That corresponds to the sudden enhancement of the SAD from the Laramie measurements and has been identified as ash-rich particles (Pueschel et al., 1994). The faster transport to the northern mid-latitudes in the models than observed may have removed most of the stratospheric particles so that the aerosol lifetime in the models is about half that observed.

In addition to different transport and microphysical mechanisms, the neglection of the Cerro Hudson eruption in August 1991 that injected about 0.75–2.0 Tg S in the form of SO₂ between 12 and 18 km (e.g., Saxena et al., 1995; Bluth et al., 1997; Neely III and Schmidt, 2016; Carn, 2022) in the simulations may partially explain the lack of the observed sulfate aerosol in the southern extratropics that we find in all model scenarios. The only exception is EMAC, which included the eruption of Cerro Hudson and nudged the meteorological variables. The importance of the Cerro Hudson eruption has therefore been evaluated with ULAQ-CCM by performing two additional simulations that consider the lower and upper estimates of the SO₂ injection in addition to the Med-22 km experiment. Significant deviations from the results of Med-22 km emerge only when including the Cerro Hudson eruption with the injection of 4 Tg SO₂ at 12–18 km altitudes (Fig. S7c, g, k–n). We observe an increase in the stratospheric sulfate burden and optical depth in the SH that better reproduces the observations for the 2 months following the Cerro Hudson eruptions. However, the shorter e-folding time of stratospheric aerosol for the extra-tropical eruption does not affect the global stratospheric lifetime and is still not sufficient to explain the lack of aerosol in the SH in the following months, which we therefore attribute to transport.

The inter-model differences may depend on numerous factors that interact with one another; this makes it hard to group models by perceived similarities, for instance a similar modal scheme, similarities in the large-scale transport, or an absence of interactive stratospheric chemistry. Laakso et al. (2022), for instance, used the same climate model (ECHAM-HAMMOZ) with two different aerosol microphysics schemes, one sectional and one modal. Even just this difference produced an effective radius up to 52 % greater in the sectional scheme than in the modal scheme simulation for the same amount of injected SO₂. Further, Niemeier et al. (2020) showed that, in two models with a similar modal scheme but different vertical advection (CESM-WACCM-110L and MAECHAM-HAM), the resulting vertical distribution of the aerosol cloud can be substantially different. Even in the same model (CESM1-WACCM), Richter et al. (2017) showed that the presence or not of interactive chemistry could strongly affect the local stratospheric warming, and thus the residual vertical velocity changes, due to feedback from the changing ozone. In our case, all of these differences are compounded; therefore it is hard to identify which exactly is the cause of the disagreement. Furthermore, in all the works cited above, SO₂ was injected continuously for a number of years rather than in an impulsive way, whereas in the case of a volcanic eruption, the synoptical conditions at the time of the eruption play an important role (Thomas et al., 2009; Toohey et al., 2014; Niemeier et al., 2021; Jones et al., 2016). In our case, the experimental protocol requires the consistency of the QBO with observations through the post-eruption period; nonetheless, there are smaller-scale processes and variability that are not reproducible by models with a coarse resolution that would affect the initial state of the system, such as the formation of a mesocyclone during the first day after the eruption (Chakraborty et al., 2009) or the passage of Typhoon Yunya within 75 km northeast of the eruption (Oswalt et al., 1996).

5 Conclusions

The ISA-MIP HErSEA experiment protocol was designed to investigate the differences and the consensus among a group of climate models, all with interactive stratospheric aerosol microphysics, by comparing them with measurements after the Mt. Pinatubo eruption in 1991. This is done through a well-defined experimental protocol with different sets of initial parameters for the stratospheric SO₂, in terms of both magnitude (5, 7, or 10 Tg S injected) and altitude of the SO₂ cloud (18–20, 21–23, 18–25 km, uniformly distributed). One important finding from this intercomparison is that there is now a general consensus among the models that an SO₂ emission amount at or below the lower end of the observed stratospheric SO₂ mass loading (14–23 Tg) is required to reproduce the observed sulfate aerosol loading from that time period. However, the set of injection parameters that best fits the observation changes in some models depending on the variables to be considered (aerosol optical depth, effective radius, sulfate burden, surface area density).

The main reason for the disagreement with observations is stratospheric transport, which is too fast towards the northern mid-latitudes for some models or results in stronger tropical confinement in others. The transport consequently influences the growth of sulfate aerosols and their global distribution, which in turn affects the persistence of aerosols in the stratosphere, with a feedback on the transport itself (Brühl et al., 2015; Niemeier and Schmidt, 2017; Visioni et al., 2018b). Other reasons could be related to the absence of processes such as the absence of the Cerro Hudson eruption in the southern extratropics 2 months after the Pinatubo eruption, which may partly explain the initial lack of sulfate aerosols in the Southern Hemisphere and the omission of ash injection, which would be crucial in the early days/months to better reproduce the initial evolution of the cloud. Our results highlight the need for some specific experiments that might be needed to disentangle the different components that contribute to the overall uncertainty. For instance, simulations that nudge stratospheric transport to reanalyses (as done in  Schmidt et al., 2018, in CESM(WACCM)) in multiple models could clarify the role of different microphysical schemes. Similarly, consistently turning interactive stratospheric chemistry on and off in multiple models could highlight the importance of ozone feedback (as done in Richter et al., 2017). Last but not least, dedicated tracer experiments of an idealized volcanic cloud excluding chemical microphysical processes are necessary to assess the role of stratospheric transport in the models.

Overall, considering the best set of initial parameters, differences between models and observations remain, and the inter-model differences are still large, as found before in other multi-model experiments of explosive volcanic clouds (i.e.,  Tambora in Marshall et al., 2018; Clyne et al., 2021). We also note that the observations themselves show disagreement, sometimes as high as inter-model differences, because of various issues with the saturation or sensitivity of the particular instrument. Our observations around the reliability of the measurements during the Pinatubo event highlight the future need for more observations in order to be better prepared for future explosive volcanic eruptions (Newhall et al., 2018; Marshall et al., 2022), both for understanding short- and long-term impacts and as a benchmark test for current Earth system models. In the absence of large volcanic eruptions in the early 21 century, where a wealth of observational data exist, it might therefore be also an alternative to focus on moderate eruptions in future aerosol model intercomparison studies, e.g., the Raikoke eruption in 2019.

As a first study of the inter-model differences within ISA-MIP HErSEA, we focused on the aerosol optical depth and the variables on which it depends, such as the loading and size of the sulfate aerosols. Therefore, we suggest for follow-up studies the comparison of radiative forcing and ozone changes, which immediately follow the analyses done here.

Appendix A: Analysis of model output

A1 Taylor diagrams

In Sect. 3.1 we use Taylor diagrams (Taylor, 2001) in order to summarize all the information regarding the reproducibility of the stratospheric optical depth simulated compared to satellite observations. Taylor diagrams provide a concise statistical summary of how well patterns from simulations and observations match each other in terms of their correlation (COR; azimuthal angle), their root-mean-square difference (RMSD; proportional to the distance between the observations (gray and black circles on the x axis) and experiments (colored circle)), and the ratio of their variances (SDs on x and y axis). SDs, RMSDs, and CORs are calculated for zonal values of the stratospheric AOD for two different time periods (first year and second year after the eruption). Therefore, similar SD, higher COR, and lower RMSD mean similar amplitudes of variation in terms of latitudinal distribution and time evolution.

A2 Effective radius

The effective radius is calculated as the ratio of the third and second moments of the number size distribution of the aerosol particles. This results in Eq. (A1) for models with a sectional scheme; in this case, the sum is over the bins, n_i is the number of particles, and r_i is the radius of particles in each bin. In models with a modal scheme, the effective radius is calculated as the sum over the modes as in Eq. (A2), where SAD_i is the surface area density, and vol_i is the volume density. In EMAC (modal scheme) the quantity is estimated from the median radii for accumulation- and coarse-mode particles since it was not stored in the output.

$$\begin{array}{}\text{(A1)}& {\displaystyle }\mathrm{reff}={\displaystyle \frac{{\sum }_i{n}_i\cdot r_i^{\mathrm{3}}}{{\sum }_i{n}_i\cdot r_i^{\mathrm{2}}}}\text{(A2)}& {\displaystyle }\mathrm{reff}={\displaystyle \frac{\mathrm{3}\cdot {\sum }_i{\mathrm{vol}}_i}{{\sum }_i{\mathrm{SAD}}_i}}\end{array}$$

The stratospheric effective radius (reff_strat) for the models and SAGE II is calculated in Eq. (A3) by integrating the provided effective radius (reff) from the tropopause to the top of the atmosphere weighted with the SAD. The thickness of the vertical layer (h) is calculated from the hypsometric equation (Eq. A4).

$$\begin{array}{}\text{(A3)}& {\displaystyle }{\mathrm{reff}}_{\mathrm{strat}}={\displaystyle \frac{{\sum }_z{\left(\mathrm{SAD}\cdot h\cdot \mathrm{reff}\right)}_z}{{\sum }_z{\left(\mathrm{SAD}\cdot h\right)}_z}}\text{(A4)}& {\displaystyle }h={\displaystyle \frac{R\cdot T}{g}}\cdot \ln {\displaystyle \frac{P_{z+\mathrm{1}}}{P_z}}\end{array}$$

For the OPC measurements, we calculate the stratospheric effective radius (Eq. A5) as in Kleinschmitt et al. (2017) for the updated UWv2.0 data set. The measurement error bars consider a 40 % uncertainty in SAD and vol and assume a correlation coefficient of 0.5 between SAD at different altitudes, vol at different altitudes, and SAD and vol at the same altitude.

$$\begin{array}{}\text{(A5)}& {\mathrm{reff}}_{\mathrm{strat}}={\displaystyle \frac{\mathrm{3}\cdot {\sum }_z{\mathrm{vol}}_z}{{\sum }_z{\mathrm{SAD}}_z}}\end{array}$$