Assessing the stratospheric temperature response to volcanic sulfate injections by Mt. Pinatubo: insights from the Interactive Stratospheric Aerosol Model Intercomparison Project

Some major volcanic eruptions, such as the one of Mt. Pinatubo in 1991, can inject large amounts of sulfur dioxide (SO₂) into the stratosphere, leading to a volcanic aerosol cloud. This dense aerosol cloud induces a radiative heating of the stratosphere, causing ozone and water vapour changes, thereby altering middle atmospheric dynamics and chemistry. The scale of these impacts for varying injection amounts and heights on stratospheric temperature anomalies is still highly uncertain. Here we analyse specially designed chemistry-climate model experiments following the Historical Eruptions SO₂ Emission Assessment Protocol (HErSEA) under the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA-MIP). The results confirm our general understanding of the stratospheric aerosol forcing due to extra SO₂ injection, while simultaneously highlighting structural differences between models. Overall, for the Pinatubo-like experiments the multi-model mean temperature anomalies agree well with meteorological reanalysis data sets, and we find that in most cases, differences between models are larger than differences for individual models across experiments with varying injection amounts and altitudes. Differences in transport, radiative transfer, and microphysics as well as the characterisation of aerosol size distributions, play a crucial role in the emergence of the spread in the modelled temperature response. Our results show further that the sensitivity of the stratospheric temperature response to model selection is also apparent in other MIPs. Hence, we argue for caution in attribution studies and the interpretation of stratospheric aerosol injection experiments relying on individual or few models.

1 Introduction

Large volcanic eruptions can affect Earth's radiation budget and climate through injection of sulfur dioxide (SO₂) and other gases into the stratosphere. Today, it is widely understood that the injected SO₂ is converted to sulfate (SO₄) aerosol within a few days to weeks, depending on the background chemical conditions, thereby thickening the stratospheric aerosol layer. This enhancement of the stratospheric aerosol layer scatters incoming shortwave (SW) radiation, absorbs solar near-infrared (near-IR) and terrestrial longwave (LW) radiation, and affects the chemical composition of the stratosphere (Robock, 2000). Due to the reduced insolation, Earth's surface is cooled, while the stratosphere is warmed through enhanced absorption by the aerosol particles formed (McCormick et al., 1995; Stenchikov et al., 1998; Lacis et al., 1992). These changes in atmospheric composition and thermal structure can have direct and indirect effects on stratospheric circulation, dynamics, and chemistry (Toohey et al., 2014; Diallo et al., 2017; Barnes et al., 2016; Zanchettin et al., 2022; Driscoll et al., 2012).

In the light of the well-documented cooling following volcanic sulfate emissions, stratospheric aerosol injection (SAI) has been discussed among a group of climate intervention strategies that might help limit global warming (Crutzen, 2006). In recent years substantial efforts have been directed to explore the feasibility of such intervention utilizing experiments with chemistry-climate models (CCMs) and Earth System Models (ESMs) (e.g., Tilmes et al., 2022; Richter et al., 2017) and analysing potential impacts on atmospheric dynamics (Aquila et al., 2014; Wunderlin et al., 2024), the environment and climate (Huynh and McNeill, 2024; Haywood et al., 2025). Although there are important differences between volcanic emissions and continuous SAI (Robock et al., 2013; Laakso et al., 2022), large volcanic eruptions remain the main natural analogue to evaluate how state-of-the-art aerosol microphysical models simulate SAI.

One of the last major volcanic eruptions that significantly affected global climate was the eruption of Mt. Pinatubo, a tropical volcano in the Philippines (15° N, 120° E), on 15 June 1991 (Robock, 2000). It was well observed by satellites (e.g., Long and Stowe, 1994; McCormick, 1987; Randel et al., 1995), airborne (e.g., McCormick et al., 1995; Pueschel et al., 1994; Borrmann et al., 1995; Deshler et al., 2003) and ground-based measurements (e.g., Nagai et al., 2010; Dutton and Christy, 1992; Antuña et al., 2002; Good and Pyle, 2004). The eruption injected SO₂, water vapour and ash into the stratosphere. There are, however, substantial uncertainties associated with injection magnitudes and height (see detailed discussion in Quaglia et al., 2023). The global mean surface cooling following the Mt. Pinatubo eruption was estimated as up to −0.5 K, and a duration of the negative temperature anomaly of 18–36 months (e.g McCormick et al., 1995; Parker et al., 1996). Some uncertainties are, however, associated with these estimates, e.g. a recent study reports a weaker cooling of −0.2 K and also a shorter duration (Boretti, 2024). In contrast to surface cooling, pronounced stratospheric warming has been reported in the literature, with reanalyses (e.g., Fujiwara et al., 2015) and observations (e.g., Randel et al., 2016; Robock, 2000) showing the largest warming signals of up to 2–3 K over a year in the tropical lower stratosphere.

Recent research (DallaSanta et al., 2019) has shown that the warming of the tropical stratosphere following a tropical (Pinatubo-like) volcanic eruption leads to a strengthening of the polar vortices and a poleward shift of the jet in both hemispheres during winter. In addition, it was shown that precipitation is indirectly altered, namely through effects of stratospheric warming on atmospheric circulation (McGraw and Polvani, 2024; Simpson et al., 2019). In terms of chemistry in the middle atmosphere, ozone concentrations have been shown to be affected both by heterogeneous chemistry on the sulfate aerosol surfaces and by dynamical effects resulting from additional heating of the stratosphere (Aquila et al., 2013; Muthers et al., 2015; Dhomse et al., 2015). Ozone concentrations can increase or decrease depending on altitude and thermal environment as SO₄ load effects heterogeneous reactions involving NO_x, HO_x, BrO_x and ClO_x and thereby alter atmospheric SW heating rates (Richter et al., 2017). Furthermore, increases in the stratospheric water vapour mixing ratio can be found due to the heating by volcanic aerosols leading to cold-point tropopause warming (Keeble et al., 2021; Kroll et al., 2021).

Over the last few decades, multi-model intercomparison projects (MIPs) have emerged as a widely used tool within the climate research community. These MIPs usually serve two purposes: (i) to provide outputs for process or trend studies under standardized protocols specifying boundary conditions and forcings, and (ii) to identify and explore biases between models and observations and inter-model differences. The Coupled Model Intercomparison Project (CMIP) supporting the Assessment Report Cycles of the Intergovernmental Panel on Climate Change (IPCC) and the Chemistry-Climate Model Initiative (CCMI and preceding MIPs) supporting the WMO/UNEP Ozone Assessment Reports constitute two of the most prominent MIPs. As these MIPs have not been designed to explicitly study the atmospheric/climatic impacts of volcanic eruptions or SAI, other MIPs such as the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA-MIP; Timmreck et al., 2018), the Model Intercomparison Project on the climatic response to Volcanic forcing (VolMIP; Zanchettin et al., 2016) or the Geoengineering Model Intercomparison Project (GeoMIP; Visioni et al., 2023) have been initiated to address these concerns.

Within the ISA-MIP activity for interactive stratospheric aerosol CCMs (Timmreck et al., 2018), an intercomparison of simulations of the Pinatubo aerosol cloud (Quaglia et al., 2023) was made, to identify any common model-observation biases. Consistent also with findings from the 2006 Assessment of Stratospheric Aerosol Properties (ASAP) report modelling chapter (SPARC, 2022) one common model bias sees models tending to over-predict the optical depth of the Pinatubo cloud, compared to the benchmark Global Space-based Stratospheric Aerosol Climatology (GloSSAC)/Stratospheric Aerosol and Gas Experiment II (SAGE-II) gap-filled dataset. Dhomse et al. (2020) suggest the discrepancy may be due to the models not including volcanic ash, the effect shown to remove 43 % of emitted sulfur in the 2014 Kelud eruption (Zhu et al., 2020). Whilst there were substantial differences among model predictions, another persistent difference between the models and observations was that the simulations removed volcanic aerosol too quickly during the decay phase after peak optical depth. This result is in agreement with the VolMIP Tambora interactive stratospheric aerosol comparison (Marshall et al., 2018; Clyne et al., 2021), which found that key processes are missing in some models, leading to large inter-model spread in predictions. Zanchettin et al. (2022) investigated the surface temperature response in VolMIP and found good agreement among the models regarding the emergent surface cooling but differences in the radiative fluxes.

To date, however, no study has systematically addressed the stratospheric temperature response to volcanic aerosol in MIPs. Here we aim to close this gap by studying the stratospheric temperature response to volcanic aerosol forming after the Mt. Pinatubo eruption in the ISA-MIP models, and illustrating and explaining important model similarities and differences across models and between models and observations.

Our manuscript is structured as follows: The experimental design, information about the participating models in the relevant ISA-MIP experiments, data sources for the additional MIPs considered, as well as reanalyses and observational data sets are provided in Sect. 2. Section 3 shows model results of SO₄ mass concentration and the stratospheric temperature response to different injection strategies, which are summarized and discussed in Sect. 4, which also presents ideas for directions of future community research.

2 Data and methods

In the subsequent subsections, we provide information on the experimental design of the ISA-MIP HErSEA experiments, building the foundation of this study, participating models, reanalysis products and observational data used to evaluate model outputs, data processing and data from further MIPs. For brevity we provide only information of direct relevance for the present study, for further details on the ISA-MIP HErSEA intercomparison we refer the interested reader to Timmreck et al. (2018) and Quaglia et al. (2023).

2.1 Model intercomparison project and experimental design

All models participating in the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA-MIP) are using interactive aerosol microphysical modules and some also interactive chemistry (see Table 1). In the framework of the Historical Eruptions SO₂ Emission Assessment (HErSEA) experiments, modelling groups performed transient Atmospheric Model Intercomparison Project (AMIP)-type runs of the Mt. Pinatubo eruption in which sea surface temperatures (SSTs) and sea ice extent are prescribed as monthly climatologies between 1991 and 1995 from the Met Office Hadley Center Observational data set (Rayner et al., 2003). Within HErSEA, five different emission scenarios were designed taking different amounts and altitudes of the volcanic SO₂ injection into account. Three of these scenarios are described by injections at a medium altitude (between 21–23 km) but different amounts of SO₂: low amount of 5 TgS (Low-22km), a medium 7 TgS (Med-22km) and high 10 TgS (High-22km). The medium-amount injection of 5 TgS has two additional scenarios with different injection altitudes: a shallow one at lower altitudes (18–20 km, Med-19km) and one over a deep altitude range (18–25 km, Med-18–25km). For brevity we refer to these experiments and the performing models hereinafter as ISA-MIP. We note in passing that EMAC model outputs are only presented in the analysis when discussing model similarities and differences in SO₄ concentrations and transport, as EMAC did not fully adhere to the ISA-MIP experimental protocol (Quaglia et al., 2023).

The ISA-MIP simulation period spans from January 1991 to December 1995 (for UM-UKCA only until December 1993) and the eruption of Mt. Pinatubo is timed on 15 June 1991. SO₂ injection is specified via a uniform vertical profile in a single grid cell at the respective altitudes close to the Mt. Pinatubo location (15° N, 120° E). The precision is defined by the specific vertical and horizontal model resolutions. The injection strategy differs only for EMAC (for which reason we flag it as EMAC^* in the figures and tables), where volcanic SO₂ is injected at a 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, Brühl et al., 2021). The eruptions of Cerro Hudson (10 August 1991), Spurr (27 June 1992), and Lascar (18 April 1993) are included in EMAC and WACCM6-CARMA but they are not included in the remaining ISA-MIP simulations.

Following the experimental protocol, all models matched the observed post-Pinatubo QBO transition, either by nudging or from re-initialising from a prior control run. This ensures a consistent QBO state across different models, minimizing the effect of model-specific QBO dynamics on model spread on aerosols. Since the main focus of the HErSEA experiment was to improve our understanding of the long-term evolution of the Pinatubo volcanic cloud (Quaglia et al., 2023; Timmreck et al., 2018), the modelling teams agreed that radiative effects of ash, water vapour and SO₂, as mainly relevant in the initial weeks post eruption (Stenchikov et al., 2021; Abdelkader et al., 2023), could be excluded. For completeness, it shall be noted that these effects can be important. For example, fine ash, due to the additional radiative heating below the sulfate cloud, can affect the initial transport direction of the aerosol cloud and inclusion of this feature can help represent modeled transport better (Niemeier et al., 2009, 2021).

For each different injection setting a three-member ensemble was submitted, except for ULAQ-CCM and EMAC, which submitted only one realization. In the respective model description sections the generation of the ensemble members is described. Unless otherwise specified, the ensemble mean of models is shown in the results.

2.2 Participating models

2.2.1 ECHAM5-HAM

MAECHAM5-HAM (further ECHAM5-HAM) is based on the ECHAM5 global climate model (GCM) (Giorgetta et al., 2006) and is used as a high-top model in the middle atmosphere (MA) version. It is also interactively coupled to the aerosol microphysical model HAM (Stier et al., 2005). Simulations were performed with a horizontal resolution of approximately 2.8°  in longitude and latitude, a spectral truncation at wave number 42 (T42), and 90 vertical layers up to 0.01 hPa (approximately 80 km). This allows for an interactive QBO simulation. The mass conservation during advection is ensured by the scheme used in ECHAM5 for transport calculations, as outlined in the work by Lin and Rood (1996). Present-day conditions are prescribed with a land-sea mask, SST, and sea ice. These values are set to climatological averages from 1950 to 2000, as described in Hurrell et al. (2008).

The HAM aerosol microphysical model (Stier et al., 2005) calculates sulfur oxidation and sulfate aerosol formation, including nucleation, accumulation, condensation, and coagulation processes. The width of the HAM modes has been adapted to conditions with a high sulfur load after volcanic eruptions. The aerosols are prescribed in three modes, nucleation, Aitken, and accumulation, with a fixed width (Niemeier et al., 2009).

HAM was also adapted to stratospheric conditions by applying simple stratospheric sulfur chemistry above the tropopause (Timmreck, 2001; Hommel et al., 2011). HAM prescribes monthly oxidant fields of nitrogen oxides, hydroxyl radicals (OH), and ozone, as well as photolysis rates of OCS, H₂SO₄, SO₂, SO₃, and O₃. Sulfate has been coupled with the radiation scheme of ECHAM for SW and LW radiation. Nucleation processes are simulated with a parametrisation after Vehkamäki et al. (2002), further adapted to high stratospheric sulfur loads by including collision processes.

Further details are described in Niemeier et al. (2021). We note, that the ensemble members were produced by increasing stratospheric horizontal diffusion by perturbing the respective parameter between 1.0, 1.0001, and 1.001 and setting it at the beginning of the simulation in January 1991.

2.2.2 ECHAM6-SALSA

The ECHAM6-SALSA refers to the aerosol–climate model ECHAM6.3-HAM2.3-MOZ1.0 (Kokkola et al., 2018) with the sectional aerosol module SALSA. The host component is the ECHAM6.3 general circulation model (Stevens et al., 2013), which is interactively coupled to the HAM aerosol module. A T63L95 resolution was used, corresponding to a horizontal grid of approximately 1.9° × 1.9°. The vertical atmosphere was divided into 95 layers, extending up to 80 km. This resolution enables the model to capture the QBO in the tropical stratosphere, which plays a key role in stratospheric aerosol transport (Laakso et al., 2022). Similarly to ECHAM5, ECHAM6 employs the flux-form semi-Lagrangian transport scheme of Lin and Rood (1996) to simulate tracer advection.

The HAM aerosol microphysical model simulates the emissions, removal processes, and radiative properties of the major aerosol species: sulfate, organic carbon, black carbon, sea salt, and mineral dust. It also includes gas- and liquid-phase sulfur chemistry. Concentrations of OH and ozone are prescribed from a monthly mean climatology. Anthropogenic aerosol emissions follow the Community Emission Data System (CEDS) for the CMIP6 emission inventory, while sea salt and dust emissions are calculated online.

In the ECHAM6-SALSA, aerosol microphysical processes – including nucleation, condensation, coagulation, and hydration are represented by the sectional aerosol model SALSA. Aerosols are resolved into 10 size bins, with the seven largest bins further separated into soluble and insoluble fractions. A comprehensive description of SALSA is provided in Kokkola et al. (2018).

Ensemble members were generated by introducing slight perturbations (0.0001 %) to one of the model tuning parameters, namely the snow formation rate by aggregation, applied only in January 1991 for each ensemble member.

2.2.3 MIROC-CHASER

The global chemistry–climate model MIROC–CHASER (Sudo et al., 2002; Watanabe et al., 2011) consists of the Model for Interdisciplinary Research on Climate (MIROC) version 4.0 and the atmospheric chemistry model CHASER (Sudo et al., 2002) and the Spectral Radiation–Transport Model for Aerosol Species (SPRINTARS) (Watanabe et al., 2011). For this study, the model is set up with a 2.8° × 2.8° horizontal resolution and 57 vertical levels up to 52 km. With this vertical resolution, QBO is not internally generated. The simulated tropical zonal winds in the stratosphere (15° N–15° S, above 90 hPa) were nudged to the observations compiled by Free University of Berlin (Naujokat, 1986) with a relaxation time of 20 d. Tracer advection was simulated using the piecewise parabolic method (Colella and Woodward, 1984) and a flux-form semi-Lagrangian scheme (Lin and Rood, 1996).

The aerosol module SPRINTARS tracks sulfate aerosol with three modes (nucleation, Aitken, and accumulation) and uses the bulk approach for black carbon and organic matter, dust, and sea salt (Sekiya et al., 2016). Nucleation is based on Vehkamäki et al. (2002). Condensation and evaporation processes are calculated with the equilibrium gaseous sulfuric acid as a function of temperature using a model referred from Kulmala and Laaksonen (1990), while the coagulation process follows the same scheme as ECHAM5–HAM (Stier et al., 2005) with van der Waals interaction as formulated by Alam (1987).

The chemical scheme in CHASER includes 93 species, as well as 263 reactions. Sulfur chemistry is included in the form of 12 reactions, as well as 4 of the main sulfur species (SO₂, SO₄, DMS, OCS) (Sekiya et al., 2016).

The initial conditions for the three ensemble members were taken from a 20 year timeslice simulation.

2.2.4 SOCOL-AERv2

Solar Climate Ozone Links version 3 with the interactively coupled aerosol model AER version 2 (SOCOL-AERv2) is an aerosol-chemistry-climate built with MA-ECHAM5 as the general circulation model. It is run in T42 configuration with a horizontal resolution of roughly 2.8°, 39 height levels going up to ∼ 80 km and a time step of 15 min. With this vertical resolution, the QBO is not internally generated, but nudged instead. The transport scheme, described in detail by Lin and Rood (1996), conserves the mass and preserves the shape of water vapour and chemical species during advection. A comprehensive model description is given in Stenke et al. (2013) and Feinberg et al. (2019). For radiative transfer calculations SOCOL-AERv2 employs the SW parameterization of Fouquart and Bonnel (1980) on 6 wavelength bands and the Rapid Radiative Transfer Model (RRTM) (Mlawer et al., 1997), which uses the correlated-k method, on 16 LW bands. Sulfate aerosol microphysics is treated by the sectional sub-model AER that calculates nucleation, coagulation, condensation, evaporation, and sedimentation of aerosol particles. It utilizes 40-size bins with radii from 0.39 nm to 3.2 µm (Sheng et al., 2015). Chemistry is represented in the module MEZON that calculates the interactions of 89 species (e.g. oxygen, hydrogen, nitrogen, …) using 234 reactions, including the detailed sulfur cycle (Brodowsky et al., 2024) and heterogeneous reactions on sulfuric, as well as polar stratospheric cloud (PSC) aerosol particles.

The time step for aerosol, chemical and radiative processes is 2 h, however, 20 sub-steps are calculated for aerosol microphysics, yielding an effective interval of 6 min. A three-member ensemble has been generated by modifying the carbon dioxide mixing ratio boundary condition by ±0.5 % for 1 month.

2.2.5 ULAQ-CCM

The ULAQ-CCM (University of L'Aquila Chemistry Climate Model) is a global-scale chemistry-climate coupled model with a horizontal lat × lon resolution of 5° × 6° (T21) and 126 log-pressure vertical levels (approximate pressure altitude increment of 568 m), from the surface to the mesosphere (0.04 hPa). The QBO is not internally resolved, but nudged to observed values (Morgenstern et al., 2017). A flux-form Eulerian explicit transport scheme is used in the model, with a minimum vertical and horizontal diffusivity to ensure numerical stability.

The chemistry module includes medium- and short-lived species grouped in families (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). This module includes seven sulfur species (OCS, CS₂, DMS, H₂S, SO₂, MSA, H₂SO₄), as well as the most important long-lived species relevant for stratospheric chemistry (N₂O, CH₄, CO, hydrocarbons, CFCs, HCFCs, HFCs, PFCs, halons).

The microphysical code for aerosol formation and growth includes a gas–particle conversion scheme, homogeneous and heterogeneous nucleation, coagulation, condensation, and evaporation (Pitari and Mancini, 2002; Pitari et al., 2016). It also includes heterogeneous chemical reactions on sulfuric acid aerosols and polar stratospheric cloud particles; both heterogeneous and homogeneous upper-tropospheric ice particle formation processes are also included (Visioni et al., 2018). 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.

The ULAQ radiative transfer module, operating online in the chemistry-climate coupled model, is a two-stream delta-Eddington approximation model used for chemical species photolysis rate calculation in UV-visible wavelengths and for solar heating rates and radiative forcing in UV-VIS-NIR bands (Toon et al., 1989). In addition, a companion broadband, k-distribution LW radiative module is used to compute radiative transfer and heating rates in the planetary infrared spectrum (Pitari et al., 2015). Absorption/scattering optical depths take into account Rayleigh scattering, absorption from O₃, O₂, NO₂, SO₂, H₂O, CO₂ and scattering/absorption from aerosol particles. Aerosol optical properties derived from the ULAQ-CCM sectional aerosol module are passed daily to the radiative transfer module, with appropriate wavelength-dependent values of Q-ext, g, and single scattering albedo, given the calculated size distribution of the particles.

2.2.6 UM-UKCA

The UM-UKCA simulations are performed using Global Atmosphere 4.0 configuration (Walters et al., 2014) of the UK Met Office Unified Model (UM v8.4) general circulation model, including the UK Chemistry and Aerosol chemistry–aerosol sub-model (UKCA). The GA4 atmosphere model has a horizontal resolution of 1.875° × 1.25° and 85 vertical levels (N96L85) from the surface to 85 km, gravity wave drag parametrization tuned to generate a realistic QBO, generated internally within the simulations. The transport scheme within GA4 is the semi-implicit semi-Lagrangian atmospheric dynamics, used to solve the non-hydrostatic fully compressible deep-atmosphere equations of motion (see Davies et al., 2005). Note that this differs from the ENDGAME dynamics used in the GA7 atmosphere model (Walters et al., 2019), the configuration used for UKESM1.0 (e.g., Sellar et al., 2019, 2020).

The UM-UKCA configuration here is adapted slightly from the basic GA4 described in Walters et al. (2014), with radiative effects from the interactive GLOMAP aerosol microphysics scheme, the GA4 model still based from the mass-based HadGEM2-ES aerosol scheme (e.g., Collins et al., 2011). The stratospheric ozone layer is also radiatively active in the simulations, from the interactive “CheST scheme” (“Chemistry of the Stratosphere and Troposphere”, see Archibald et al., 2020). The configuration of the GLOMAP aerosol microphysics module includes the adaptations for stratospheric aerosol described in Dhomse et al. (2014), with additionally heterogeneous nucleation of sulfate aerosol from meteoric smoke particles, as described in Brooke et al. (2017), and used equivalently in the GA4 UM-UKCA simulations in Marshall et al. (2018, 2019) and Dhomse et al. (2020).

The boundary conditions for the transient GA4 UM-UKCA Pinatubo simulations is described in Dhomse et al. (2020), using Greenhouse gas (GHG) and ozone-depleting substance (ODS) concentrations from Ref-C1 scenario from CCMI-1 (Morgenstern et al., 2017). The simulations were in atmosphere-only mode, specifying time-varying SSTs and sea ice concentrations through the post-Pinatubo period from the AMIP simulations for CMIP6 (obtained from https://esgf-node.llnl.gov/projects/cmip6/, last access: 25 March 2021).

The initialisations for the three ensemble members were selected from 20-year time-slice simulations that gave a QBO transition through 1991, approximately matching that of the ERA-Interim reanalysis (Dee et al., 2011), see Appendix A of Dhomse et al. (2020).

2.2.7 WACCM6-CARMA

In this study, the Community Earth System Model version 2.2 (CESM2.2; Danabasoglu et al., 2020) is configured with the atmospheric component of the Whole Atmosphere Community Climate Model, incorporating Middle-Atmosphere chemistry (WACCM6-MA; Gettelman et al., 2019; Davis et al., 2023), at a horizontal resolution of 1.9° × 2.5°, using the finite volume dynamical core (Lin and Rood, 1997). WACCM6-MA the high-top configuration of CESM, featuring 70 vertical levels and extending up to approximately 150 km. The QBO is not internally generated, but the tropical zonal winds are nudged to idealized observations between 4 and 86 hPa and 10° S to 10° N (Matthes et al., 2010) following Davis et al. (2023). The MA model includes comprehensive chemistry for the stratosphere, mesosphere, and lower thermosphere, while employing a simplified chemical mechanism in the troposphere. The scheme includes 59 chemical species, 217 gas-phase reactions, and 17 heterogeneous reactions.

Aerosols in both the troposphere and stratosphere (including black carbon, primary organic matter, secondary organic aerosols, sea salt, dust, and sulfate) are simulated using the sectional microphysics scheme of the Community Aerosol and Radiation Model for Atmospheres (CARMA; Yu et al., 2015; Tilmes et al., 2023). In this configuration, CARMA includes two aerosol groups: one representing internally mixed particles and another representing pure sulfate particles. Each group is resolved into 20 size bins, ranging from 0.05 to 8.7 µm for the mixed group and from 0.3 nm to 1.3 µm for the pure sulfate group. The model accounts for binary homogeneous nucleation of sulfuric acid and water, applied exclusively to the pure sulfate group using the Zhao scheme (Zhao and Turco, 1995), as well as condensational growth and evaporation of sulfuric acid for both groups. Coagulation processes are included within each group and between groups. Both aerosol groups are also subject to dry deposition, calculated using the Zhang et al. (2001) parameterization, and gravitational settling, with settling velocities determined according to Seinfeld and Pandis (1998), and wet deposition (see Tilmes et al., 2023).

Aerosol optical properties are computed using a core-shell Mie scattering model, which captures the layered structure of internally mixed particles in the mixed aerosol group, following the method of Yu et al. (2015). In this setup, black carbon and dust serve as the core materials, while the shell is composed of water and water-soluble components such as organics, sea salt, and sulfate. For the pure sulfate group, the Mie calculations assume homogeneous particles with no core and use the sulfuric acid weight percentage to determine the refractive index.

The three ensemble members were initialized from a spin-up simulation, using different years (1979, 1980, and 1981) for the initialization of the start date in the year 1980.

2.2.8 EMAC

EMAC consists of the general circulation model ECHAM5 and the Modular Earth Submodel System MESSy V2.52 (Jöckel et al., 2010) including Atmospheric Chemistry. We use a setup in the spectral resolution T63 (1.9°) with 90 layers up to 1 Pa (80 km) where the meteorology of the troposphere up to 100 hPa is nudged to the reanalysis ERA-Interim (now ERA5). SSTs and sea ice cover (SIC) are prescribed by the reanalysis. The QBO is internally generated but slightly nudged to observations compiled by the Free University of Berlin (relaxation time 2 months, based on Giorgetta and Bengtsson (1999). Species are transported with the semi-Lagrangian advection scheme of ECHAM5 (Lin and Rood, 1996). EMAC uses the radiative transfer of ECHAM5 with corrections of Thomas (2008) with 4 solar and 16 terrestrial IR bands. EMAC contains comprehensive gas-phase and heterogeneous chemistry. The applied aerosol module GMXE (Pringle et al., 2010) accounts for seven modes using log-normal size distributions (nucleation, soluble and insoluble Aitken, accumulation, and coarse modes, adjusted for stratospheric applications (Schallock et al., 2023)). Optical properties for the types sulfate, dust, organic carbon, black carbon, sea salt, and aerosol water are calculated using Mie-theory-based lookup tables for each mode consistent with the simulated size distribution of each mode. The resulting total optical depths, single-scattering albedos, and asymmetry factors are used in radiative transfer calculations, which feed back to atmospheric dynamics (especially ascent and descent).

The results from EMAC were taken from an existing 30-year transient simulation for comparison with satellite observations (Schallock et al., 2023). The 3D plume for sulfur injection by Mt. Pinatubo was derived from several weeks of SAGE-II data, leading to a spread over a latitude and altitude belt (ca. 20°, 3 km).

Table 1 Main characteristics of models participating in ISA-MIP.

2.3 Observational data and reanalysis

2.3.1 REMAP-GloSSAC-2023

To assess the model-predicted heating of the stratosphere, we first analyse how the models differ in terms of volcanic sulfate, compared to a best-estimate observational dataset (REMAP-GloSSAC-2023, https://doi.org/10.3929/ethz-b-000713396, Jörimann et al., 2025). This was derived from the observational composite GloSSAC version 2.22 (GloSSACv2.22; NASA/LARC/SD/ASDC, 2022), which provides aerosol extinction coefficients on up to four wavelength bands. From these data, a unimodal log-normal aerosol size distribution with three free parameters can be derived. This was done by iteratively checking the error between theoretical and measured extinction coefficients for the entire realistic parameter space and selecting the best-fitting parameters, following the Retrieval Method for Aerosol optical and physical Properties in the stratosphere (REMAPv1; Jörimann et al., 2025). Thus, we produced monthly stratospheric fields of H₂SO₄ number density. These fields contain useful values up to 80° latitudinally and up to 39.5 km; however, GloSSAC provides only a high-altitude climatology above 27 km, as opposed to transient data, mostly based on the well-established SAGE II and III data. Earlier versions of REMAP were also used to create stratospheric aerosol forcing data for both CMIP6 and CCMI-2022 (Sect. 2.4).

2.3.2 Reanalysis

For evaluating the temperature response in the model simulations we use data from three atmospheric reanalyses, including the European Centre for Medium-Range Weather Forecasts (ECMWF) Fifth Reanalysis of the Atmosphere (ERA5; Hersbach et al., 2020), the Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA2; Gelaro et al., 2017), the Japanese Reanalysis for Three Quarters of a Century (JRA55; Kobayashi et al., 2015). We use zonal temperature data on pressure levels (Martineau et al., 2018) from the Atmospheric Processes And their Role in Climate (APARC) Reanalysis Intercomparison Project (A-RIP) (SPARC, 2022) in this study. We are using the common grid version of this dataset, for which all diagnostics were interpolated horizontally onto a regular 2.5° × 2.5° grid for a subset of pressure levels that are common among all included reanalyses (Martineau et al., 2018).

2.3.3 Anomaly calculation

For SO₄ concentrations, absolute values can be used instead of anomalies, as the background values are very low. The temperature signal, however, is influenced by various factors, such as the QBO, seasonal fluctuations, as well as large-scale phenomena such as the El Niño–Southern Oscillation (ENSO). Therefore, anomalies are better suited to illustrate the temperature response to the Mt. Pinatubo eruption.

For reanalyses multiple linear regression (MLR) following the approach by Fujiwara et al. (2015) is used, including two QBO indices (30 and 50 hPa monthly mean zonal wind time series around the equator), the observed Niño 3.4 index lagged by 4 months for the ENSO variability, the solar 10.7 cm flux data for the solar variability and linear trend. The effect of MLR compared to anomalies with respect to the mean of the period 1980–2000 (21 years centred around the Mt. Pinatubo eruption) can be seen in Fig. A1. While the median of the individual reanalysis distributions remains unchanged, the spread decreases substantially due to the MLR.

For the models, anomalies are built using their respective Control simulations. Since the ocean is not run interactively and QBO is nudged or re-initialized at the time of the eruption, their effects should be at least partially accounted for. Models with available wind data (ULAQ-CCM, MIROC-CHASER and SOCOL-AERv2) were used to analyse the QBO contribution to the temperature anomaly signal. This was achieved by comparing the anomalies (experiment minus Control) without QBO treatment and with QBO regressed out (model_QBO in Fig. A2). We can see that the effect is much smaller compared to the results for reanalysis using the mean of the period 1980–2000 (see Fig. A1), and therefore we do not apply MLR prior to calculating the anomalies of the models hereinafter.

2.4 Further MIPs

As noted above, the REMAP method which we use to derive mass-related stratospheric aerosol properties from observations is a successor to the algorithms used in CMIP6 and CCMI-2022 forcing preparations. After retrieving stratospheric aerosol size distributions the preparation of the forcing data sets for these MIPs included calculating not only the physical aerosol properties (effective radius, number density, and surface area density) but also the optical properties: aerosol extinction coefficient, single-scattering albedo and asymmetry factor – tailored to each model's spectral resolution. For CMIP6 these data were created with a similar best-fit algorithm – using, however, an earlier version of GloSSAC (v1.1) and selecting three out of four SAGE II wavelength channels. The data set was therefore termed SAGE-3λ and is described in detail in the Supplement to Jörimann et al. (2025). The forcing for CCMI-2022 was then based on a further extended and updated version of GloSSAC (v2.0). A key difference in the REMAP code applied to the different products was a parameterization used for regions, where data was scarce.

2.4.1 CMIP6-AMIP

We use monthly temperature data of 20 Climate Model Intercomparison Project Phase 6 (CMIP6) models (Eyring et al., 2016) within the historical AMIP experiment as listed with respective ensemble sizes in Table A1. The AMIP simulations span the period 1979–2014 and use observed SSTs and SICs, as well as prescribe other external forcings including volcanic aerosols, solar variability, GHG concentrations, and anthropogenic aerosols.

2.4.2 CCMI-2022

We use monthly temperature data from the nine CCMs participating in the historical simulations (REF-D1) of the Chemistry Climate Model Initiative 2022 (CCMI-2022; Plummer et al., 2021), which are listed together with the ensemble size of each model in Table A1. Only ensemble members created through different realisations and not by perturbing initialization, physics or forcing data were selected. The REF-D1 historical hindcast simulations use forcings that follow the observed historical conditions closely and cover the period 1960–2018. The simulations have prescribed monthly mean SSTs and SIC from the HadISST data set (Rayner et al., 2003). Even though the QBO should be nudged in the CCMI-2022 ensemble, Benito-Barca et al. (2025) found that in some models QBO is not nudged but internally generated and thus, differences between observations and models can evolve from non-synchronized QBO phases. Stratospheric aerosol surface area density (SAD) is prescribed from the same REMAP methods used to relate the size distribution to the 3λ aerosol extinctions from the SAGE-II satellite measurements (see Revell et al., 2017; Jörimann et al., 2025).

3 Results

Given that the experimental design is designed to mimic the eruption of Mt. Pinatubo we start with an analysis of the time evolution of temperature anomalies in the tropics (Sect. 3.1). Thereafter, we expand the analysis to the space-time propagation of the volcanic aerosol plume (Sect. 3.2), the associated temperature response across the middle atmosphere (Sect. 3.3), and discuss the model sensitivity to changes in SO₄ concentrations (Sect. 3.4). We progress by contrasting our findings for ISA-MIP with those for other recent MIPs (Sect. 3.5) and conclude with Sect. 4.

3.1 Tropical temperature anomalies

We start by illustrating in Fig. 1a the time evolution of the mean temperature anomaly in the tropical lower stratosphere normalized to June 1991 for the five HErSEA-experiments from the ISA-MIP project along with the multi-reanalysis mean (comprising ERA5, MERRA2 and JRA55). The stratospheric heating due to volcanic aerosols, as widely discussed in the literature (e.g., Robock, 2000; Stenchikov et al., 1998; McCormick et al., 1995), can be clearly seen in the reanalysis mean. After a short ramp-up period the heating peaks in September 1991 and stays relatively constant at approximately 2 K until January 1993. Thereafter, the tropical heating signal quickly diminishes and vanishes during boreal summer (JJA) 1993.

While the overall evolution of the temperature anomaly identified in the reanalyses is broadly reproduced by the multi-model mean across all experiments, the timing and persistence of the temperature anomaly peak differ from the pattern observed. Among the individual models the temperature anomaly signal fades faster than observed (already between June and October 1992) in ECHAM6-SALSA, ULAQ-CCM and MIROC-CHASER in all experiments. Regarding the magnitude of the temperature anomaly, we find a good agreement between reanalysis and the multi-model mean for all the experiments except the High-22km. In this experiment, the injection of 10 Tg of sulfur in the form of SO₂ results in an overestimation of the temperature signal in most of the models and thus in the multi-model mean. For the other experiments we find, despite good agreement between the reanalyses and the multi-model mean, large differences between individual models. The envelope in temperature anomalies ranges between ECHAM6-SALSA, which shows the smallest temperature anomaly and the least difference across experiments and SOCOL-AerV2, which shows the strongest temperature signal and pronounced sensitivity to different injection magnitudes and altitudes.

Given the similarities at the 50 hPa level in Low-22km, Med-22km, Med-19km and Med-18-25km we expand in Fig. 1b the analysis to the vertical temperature anomaly in the tropical region for the time period with the strongest signal across models, which corresponds also to the initial peak in the reanalysis mean (September to December 1991). During this period, the reanalysis shows a peak in temperature anomaly of around 2 K between 30 and 40 hPa. While the multi-model mean captures the altitude of the peak and overall vertical distribution of the temperature anomalies well in most experiments, the magnitude is overestimated in all but Low-22km. This overestimation in the multi-model mean is due to the majority of models showing larger anomalies than observed, with several models exceeding the 3σ confidence interval of the reanalysis. With regard to the altitude of the peak temperature anomaly, experiment Med-19km marks an exception. In this experiment the anomaly peaks already around 50 hPa which can be explained by the low injection height. In summary, Fig. 1 shows that the multi-model mean, for most experiments, agrees well with the reanalysis mean in both time evolution and vertical distribution. Additionally, we find that differences between the models are larger than differences of individual models across experiments.

Therefore, to focus our study, we selected one of the five experiments as the core experiment for further analysis hereinafter, and retained the additional experiments for the discussion of important differences if they emerge. To select this experiment we used the Root Mean Square Error (RMSE) between the multi-model mean and reanalysis for the time evolution analysing the pressure levels with the strongest signal (see Table 2). As High-22km is overestimating the observed temperature response and the signal fades too fast in Med-18-25km, both experiments yield high RMSEs. While an overall best-agreeing experiment is not clearly emerging, the lowest mean RMSE is found for Low-22km and Med-22km. Given that Med-22km agrees better with reanalysis for most of the models than Low-22km (see Table A2), we selected Med-22km as our core experiment.

Figure 1 Time evolution of the observed (multi-reanalysis mean) and model-derived tropical mean temperature anomaly normalized to June 1991 (a) and tropical vertical temperature anomaly for the period September to December 1991 (b). In all panels, the multi-model mean is provided as a blue bold line, and the grey shading indicates the ±3σ confidence bound for the reanalysis mean (black, bold).

Table 2 RMSE for the time evolution of the observed (multi-reanalysis mean) and model-derived tropical mean temperature anomaly normalized to June 1991 for different pressure levels.

3.2 Space-time evolution of SO₄ concentrations

Having set the scene by illustrating the large differences in temperature signal of the Mt. Pinatubo eruption across models, we turn to the model representations of the associated forcing. As discussed above, the stratospheric heating results mainly from the absorption of terrestrial LW and solar near-IR radiation. If the absorbing particle is much smaller than the wavelength, the absorption of terrestrial LW radiation and its stratospheric heating is proportional to the mass concentration of the sulfate aerosols (Van de Hulst, 1981; Murphy et al., 2021; Lacis, 2015). In the case of volcanic aerosols, the size distribution is comparable to near-IR wavelengths; thus, differences between modelled heating can originate from both mass concentration and size distribution and need to be considered in the subsequent analysis, even though the latter has less impact on the total absorption (Ramachandran et al., 2000; Stenchikov et al., 1998, 2021). If the aerosol radius is larger than ∼ 0.5 µm the absorption efficiency starts to increase strongly with increasing aerosol size (Laakso et al., 2022).

Quaglia et al. (2023) analysed the evolution of stratospheric aerosol optical depth, stratospheric sulfate burden, and aerosol properties (effective radius and SAD) in the ISA-MIP models. The main reason for disagreement among the models analysed in that study was differences in stratospheric transport. ECHAM6-SALSA, ECHAM5-HAM, and SOCOL-AERv2 showed too fast transport from the tropics to high northern latitudes, while for UM-UKCA, transport was confined to the NH. Transport in ULAQ-CCM was found to be too confined to the tropics, likely affected by a different vertical advection scheme in comparison to the other models based on the same dynamical core, ECHAM5 or ECHAM6. Additional factors influencing the tropical confinement were differences in vertical and horizontal resolutions and the definition of the tropical pipe. Overall, most of the models showed a preference for transport to the NH. Exceptions were EMAC, which has been nudged to reanalysis to represent observed tropospheric meteorological conditions, and UM-UKCA when emissions have been spread between 15° N and the Equator instead of point injection. These models showed a better agreement with observations regarding transport to the SH. For the other models, Quaglia et al. (2023) state as reasons for the lack of aerosol (transport) in the SH, the omission of (i) ash injection, which would be crucial in the early days/months to better reproduce the initial evolution of the aerosol cloud, as well as, (ii) aerosol emissions emerging from the Cerro Hudson eruption in August 1991 in the southern extratropics. Below we build on the findings of Quaglia et al. (2023) and extend their analysis with a focus on the time evolution and vertical propagation of SO₄ concentrations and resulting temperature anomalies.

In Fig. 2 we show the time evolution of the zonal mean of SO₄ concentrations exceeding 1 µg m⁻³ across the upper atmosphere. For this comparison, the REMAP-GloSSAC-2023 data serves as an observational reference. REMAP-GloSSAC-2023 shows increasing SO₄ concentrations post-eruption in both hemispheres, which peak in the tropics until autumn 1991, and fade over the following year. In the models, zonal mean SO₄ concentrations peak earlier (in late summer/early autumn 1991) than observed. Furthermore, the initial distribution differs widely among models. While for ECHAM5-HAM and ULAQ-CCM the SO₄ concentrations are – similarly as for REMAP-GloSSAC-2023 – confined within a small and narrow fraction of the upper atmosphere, the SO₄ plume in SOCOL-AERv2 and WACCM6-CARMA spreads across large parts of the lower and middle stratosphere. These initial differences across models can be attributed to the combined effect of SO₂ conversion rates, transport and model vertical velocity (Niemeier et al., 2020). While some models (EMAC, ECHAM6-SALSA and SOCOL-AERv2) already start transporting SO₄ through the shallow branch of the Brewer-Dobson Circulation (BDC) into higher latitudes shortly after the eruption, others, such as ECHAM5-HAM, ULAQ-CCM, and UM-UKCA, have SO₄ more confined in the tropical pipe during the first months. By November 1991, all models show transport into higher latitudes; however, there are stark differences in transport pathways. Some models (e.g. EMAC, ECHAM6-SALSA, SOCOL-AERv2) show preferred transport through the shallow branch while others (e.g. UM-UKCA, WACCM6-CARMA, MIROC-CHASER) mainly transport through the deep branch until winter 1991/92. A special case is MIROC-CHASER for which SO₄ concentrations are quickly transported into the northern subtropics (20–30° N) and then further towards the pole, with limited sedimentation and no apparent subtropical transport barrier.

Despite these model differences several similarities emerge. As discussed above, a preferred meridional transport to the NH can be identified in all models but EMAC, and UM-UKCA if emissions are spread between 15° N and the Equator (see Fig. A3), in agreement with Quaglia et al. (2023). For all ISA-MIP model realizations a net descent of SO₄ concentrations begins in November 1991, likely reflecting the increased effective radius, and downward propagation of the SO₄ plume emerges between late winter and early spring 1992. From May 1992 onward, the peak in SO₄ concentration can be found over the northern polar region in all models, and only low SO₄ values remain in the lower latitudes. While this is in broad agreement with REMAP-GloSSAC-2023 observational composite retrievals, the latter shows stronger transport to the southern hemisphere than the models. Despite the northern hemispheric transport preference of the models, both, low resolution models such as SOCOL-AERv2 and intermediate resolution models such as ECHAM6-SALSA and EMAC, perform well in comparison to REMAP-GloSSAC-2023 observations regarding the magnitude and transport patterns of SO₄ concentrations. However, for ECHAM6-SALSA, the aerosol distribution between the hemispheres varies considerably across individual ensemble members, a feature highlighted in previous work (Quaglia et al., 2023; Laakso et al., 2016) pointing to strong sensitivities to the dynamic state at the injection grid point at the time of SO₂ release.

Most importantly, Fig. 2 makes large differences in aerosol transport and concentrations across models apparent; thus, we investigate next how these relate to temperature signals and differences therein.

Figure 2 Vertical evolution of zonal mean SO₄ concentrations in experiment Med-22km for individual models, from July 1991 to August 1992 (bi-monthly; columns) in comparison with REMAP-GloSSAC-2023 observations (first row). Note, EMAC is flagged with ^* as it did not fully adhere to the ISA-MIP protocol, see Sect. 2.2.8.

3.3 Comparison of SO₄ concentrations and temperature anomalies

The largest post-Pinatubo SO₄ concentrations occurred in the tropics (e.g., McCormick and Veiga, 1992) and we first assess the progression of each model's volcanic sulfate in the 20° S–20° N region. To allow for a direct comparison of SO₄ forcing and temperature anomaly we provide model-specific SO₄ concentration contour lines in all panels of Fig. 3. First, we note that for both, the multi-reanalysis mean and REMAP-GloSSAC-2023 observations, a strong relationship between the progression of peak temperature anomaly and SO₄ concentration appears in the spatial and temporal analysis. Second, in observations and in all models, the strongest temperature response occurs shortly after the eruption simultaneously with the highest concentrations in SO₄. Third, the temperature signal decays over time as the SO₄ plume propagates. Fourth, the strong correlation between temperature anomaly and SO₄ concentration identified in observations is also present in all models; however, the magnitude, timing, temporal and vertical extent of the signal differ substantially between them.

For one thing, MIROC-CHASER and ECHAM6-SALSA stand apart from the other models and the observations: MIROC-CHASER as it does not show a descending temperature signal, and ECHAM6-SALSA as it does not show any significant temperature anomalies in Med-22km. The absence of a significant anomaly in ECHAM6-SALSA might be explained by a strongly varying aerosol distribution across the individual ensemble members (see discussion above). We note in passing that in experiments resulting in a stronger temperature signal, some significant anomalies emerge (see Fig. A4). Also, ULAQ-CCM differs markedly from the remaining models given the much faster descent and fade out of the temperature anomaly while its SO₄ anomaly remains. From the remaining four models, the UM-UKCA ensemble mean appears to resemble the magnitude and pattern of temperature anomaly identified in the reanalysis. ECHAM5-HAM, SOCOL-AERv2 and WACCM6-CARMA agree with the observations regarding pronounced temperature anomalies but show in parts marked differences in magnitude, propagation and vertical extent. We attribute the deficiencies in reproducing observed patterns not solely to model issues in transport. This is exemplified as follows: while ECHAM6-SALSA and SOCOL-AERv2 agree well with reanalysis regarding the transport of the SO₄ concentration plume (see Fig. 2), they over- (SOCOL-AERv2) or underestimate (ECHAM6-SALSA) the corresponding temperature anomalies.

Figure 3 Time evolution of the vertical temperature anomalies in the tropics (20° S–20° N) in experiment Med-22km for individual models. In all panels: SO₄ concentrations are shown as contours (levels are 0.5, 1, 2, 5 µg m⁻³), dots indicate non-significant areas (p-value >0.05) for reanalysis and models with the exception of ULAQ-CCM where only one ensemble member is available, the thick black vertical line represents the Mt. Pinatubo eruption time while the thin vertical lines embed the time period September to December 1991. REMAP-GloSSAC-2023 SO₄ concentrations are shown as contours in the reanalysis panel.

To investigate inter-model differences in the simulated Pinatubo aerosol vertical profile and associated temperature effects, we focus on the time period of the strongest temperature signal-to-noise ratio across models (September to December 1991) in Fig. 4. Starting with the observations, reanalyses highlight that the peak SO₄ concentrations and temperature anomalies reside in the tropics (20° S–20° N) and across the lower and middle stratosphere (70–20 hPa). Further, the propagation of the volcanic signal over time can be clearly seen, i.e. hemispheric transport and descent (compare with Fig. 3).

Figure 4 Vertical zonal temperature anomaly in experiment Med-22km for the period September to December 1991 for individual models. In all panels: SO₄ concentrations are shown as contours (levels are 0.5, 1, 2, 5 µg m⁻³), dots indicate non-significant areas (p-value >0.05), and the vertical lines show embed the tropical region (20° S–20° N). REMAP-GloSSAC-2023 SO₄ concentrations are shown as contours in the reanalysis panel.

When probing the models, the observed relationship of the temperature anomalies and SO₄ concentrations is present in all. Somewhat an exception is ULAQ-CCM, for which the location of the largest temperature anomalies can be found below the maximum of the SO₄ concentrations. One may attribute this to a larger effective radius (exceeding 0.5 µm in radius in the tropics) in ULAQ-CCM compared to the other models in the lower stratosphere. However, given that WACCM6-CARMA also shows an tropical effective radius exceeding 0.5 µm while still smaller than ULAQ-CCM (see Fig. A6), such interpretation may not hold also pointing to limitations of effective radius in explaining the location of temperature anomalies. Also, MIROC-CHASER differs from the remaining models as its maximum SO₄ concentrations and thus peak temperature anomalies emerge due to too strong transport (see Fig. 2) outside of the tropics. For the remaining models, the spatial structure of the SO₄ concentrations agrees quite well with the reanalysis data. However, temperature responses across these models are very different in magnitude. Given the similarities in the SO₄ concentrations, we do not attribute these differences to transport but rather to model specifics in radiative transfer and on the simulated particle size distributions and thus consequently atmospheric heating rates.

In addition, the vertical structure and altitude of aerosol-driven temperature anomalies differ within the model ensemble. While the majority of the models show the most pronounced tropical warming signal centred around approximately 40 hPa (while expanding between 70 and 20 hPa), the peak anomaly for ULAQ-CCM is located substantially lower. The peak temperature response in WACCM6-CARMA and MIROC-CHASER emerges primarily at higher altitudes (see Fig. 4). In summary, the magnitude of the temperature anomaly corresponds to the spatial structure of SO₄ concentration for individual models, but comparison across models shows that the magnitude of the obtained warming differs widely. We note that the magnitude of warming also varies across experiments, with the largest (smallest) signal in High-22km (Low-22km), see Fig. A5. Therefore, we investigate next the model-specific warming per µg m⁻³ of SO₄, taking all the experiments into account.

3.4 Model sensitivity to SO₄ forcing

In this section we analyse the model sensitivity to SO₄ forcing. This metric enables us to assess potential structural differences among the models in the magnitude of stratospheric heating in response to differences in simulated sulfate aerosol. The normalisation to simulated sulfate enables the assessment solely of differences in the sensitivity of models, avoiding the problem of there also being differing heating effects from spatio-temporal differences in simulated SO₄ mass. Specifically, to calculate the post-Pinatubo sensitivity, we regress the temperature anomaly for each grid point over time against SO₄ concentration and average the regression coefficients with p-values smaller than 0.05 across the tropical stratosphere (see Fig. 5). Having previously addressed ensemble means, we apply this method to all ensemble members individually. Further, we use all available experiments with varying injection height and amount in this analysis to increase data availability. This allows us to build the relationship between SO₄ concentrations and resulting temperature anomalies considering the full range of responses for each model. As an observational benchmark we use REMAP-GloSSAC-2023 SO₄ concentrations to regress the temperature signal in the reanalysis data. To be consistent with all previous analyses, we use the pressure range in which the strongest signal in observations and most models is observed (70 to 20 hPa, see Fig. 4). Focusing on the median sensitivity across models, Fig. 5 reveals large differences, from hardly any sensitivity to a response of up to 2 K per µg m⁻³. It further appears that the spread in individual models broadens with increasing sensitivity to SO₄, i.e., models with low sensitivity provide a small temperature response quasi-independent of SO₄ concentration, while models with high sensitivity cover a broad range of temperature response dependent on SO₄ availability.

For example, ECHAM5-HAM (median: 1.72) and SOCOL-AerV2 (1.02) show a particularly high sensitivity, also larger than estimated from the reanalyses (0.63). This high sensitivity, however, has a different root in each model: for ECHAM5-HAM it emerges from an underestimation of the SO₄ concentrations, while in SOCOL-AERv2 a high sensitivity emerges from an overestimation of the temperature response. In contrast, ECHAM6-SALSA (0.11), ULAQ-CCM (0.22) and UM-UKCA (0.28) show lower sensitivity to SO₄, but also for different reasons. As discussed above, in ULAQ-CCM the altitudes of the peak temperature signal and peak SO₄ concentrations do not align as well as for the other models. This mismatch results in regressing a mixed signal that falls to lower sensitivity. For ECHAM6-SALSA and UM-UKCA the lower sensitivity can be explained by the weak temperature response despite the well-reproduced SO₄ concentrations. The sensitivities of MIROC-CHASER (0.88), and WACCM6-CARMA (0.70) complete the analysis and are those that lie closest to the reanalysis estimate. We note that these sensitivities need to be interpreted with caution as they do not characterize the pure temperature response to the SO₄ concentration, i.e. the model-specific estimate also includes the contributions from aerosol size distribution and the effective heating emerging from the radiative transfer scheme, and QBO effects.

Figure 5 Stratospheric temperature sensitivity to SO₄ absorptive-heating of the ISA-MIP models and reanalysis in the tropics (20° S–20° N) between 70 and 20 hPa. The boxplots contain data for altitude weighted longitude – latitude pairs, obtained by regressing temperature anomalies against SO₄ concentrations over time.

3.5 Comparison to other Multi-model Intercomparison Projects

These results show that even when following a strict simulation protocol, as put forward by ISA-MIP HErSEA, model outcomes for aerosol concentration and subsequent temperature change in the middle and lower stratosphere differ widely. Some of these model differences may be attributed to the use of specific aerosol schemes and/or interactive chemistry in some of the models (see Table 1). Thus, we expect that the inter-model spread can be reduced if even stricter protocols are applied. We test this hypothesis through recent MIPs under the auspices of the World Climate Research Programme (WCRP) that use prescribed SADs: CCMI-2022 REF-D1 and CMIP6-AMIP. While for the ISA-MIP models a Control run without the Mt. Pinatubo eruption is available and used for computing forced anomalies, no such Control is available for CCMI-2022 and CMIP6-AMIP. Thus, for these two MIPs, and also for the reanalyses in this comparison, we compute anomalies with respect to the mean of the period 1980–2000 (21 years centred around the Mt. Pinatubo eruption). We note in passing, that such anomaly calculation yields only a slightly broader range of anomalies for the reanalysis than the MLR calculation used above (see Fig. A1), which is, however, still substantially narrower than for the ISA-MIP models, see Fig. 6. We focus on the time period, latitude band and pressure range in which the strongest signal in observations is found (September to December 1991, 20° S–20° N, 70 to 20 hPa).

As in our MLR analysis, we find that the temperature anomalies derived for ERA5, MERRA2, and JRA55 agree well with each other and also show a similar spread for the selected period and vertical range (see Fig. 6). Our analysis also shows that the median/mean of the multi-model temperature anomaly agrees well with reanalyses for all MIPs considered, while the distribution of the temperature anomaly is much wider for MIPs than for reanalyses, indicating that the agreement in the median/mean occurs through error compensation when averaging across models.

For CCMI-2022 we find a slightly lower magnitude of temperature anomaly compared to ISA-MIP, but a median still in close agreement with reanalyses. More interestingly, however, the intra-ensemble spread (evaluated from minimum to maximum values) in CCMI-2022 is even increased (longer whiskers, wider box) compared to ISA-MIP. We interpret this outcome as follows: model complexity (i.e. inclusion/absence of interactive aerosol microphysical schemes) appears to be of secondary importance compared to other factors such as the radiative transfer scheme or model selection. Expansion of this analysis to CMIP6-AMIP further supports this interpretation as for these models the intra-model spread is even larger (longer whiskers, narrower box) while the median is marginally larger than for ISA-MIP models or reanalyses. The hypothesis that there is strong model-dependency of the temperature anomalies obtained is further supported when for the same comparison only models with more than one ensemble member are taken into account (see light orange, light green and light purple boxplots in Fig. 6, respectively). On the one hand, an increase or decrease in the multi-model median is found for individual MIPs when subset, on the other, also the intra-model spread is altered. For convenient reference we indicate individual models with dots per MIP (and subset), which further corroborates the sensitivity to model selection as visual inspection clearly shows that selecting or excluding a model (or subset of models) would strongly affect the overall result.

Figure 6 Temperature anomaly distributions for reanalyses (left), ISA-MIP and CCMI-2022 (center) and CMIP6-AMIP (right) in the tropics (20° S–20° N) averaged between 70 and 20 hPa in the period September to December 1991. The number in brackets indicates the number of models used. Solid vertical line separates reanalysis from models. Light orange, light green and light purple boxplots represent models with more than one ensemble member. Individual model results are provided as dots in the same colour of the corresponding MIP/sample. For convenient reference we note the order of ISA-MIP models in descending temperature anomaly: SOCOL-AERv2, ECHAM5-HAM, WACCM6-CARMA, ULAQ-CCM, UM-UKCA, ECHAM6-SALSA, MIROC-CHASER.

4 Discussion and conclusions

Here we use the well-defined experimental protocol of the ISA-MIP HErSEA experiments to investigate the stratospheric temperature response to the Mt. Pinatubo eruption of 1991. To this end, we compare outputs of a group of climate models with interactive stratospheric aerosol microphysics to observations and reanalysis. We explicitly focus on the evolution of the temperature response to changes in SO₄ concentration to highlight similarities and differences across the models under uniform conditions for volcanic SO₂ injection in magnitude and altitude. We find that the multi-model mean of stratospheric temperature anomalies agrees well with the reanalysis mean for these experiments (±1 K in the first four months after the eruption, see Fig. 1). Individual ISA-MIP HErSEA experiments differ in injection amount and altitude. Overall, the model results confirm our understanding of stratospheric aerosol forcing, i.e. that larger injection of SO₂ leads to a stronger warming response. Injection at a lower altitude leads to faster removal of the volcanic aerosols from the stratosphere and therefore a shorter duration of the temperature signal compared to reanalysis. While these physical relationships are present in most of the models, substantial spread in magnitude and duration of the temperature response is apparent, a result somewhat surprising given the strict experimental protocol. We find that differences between the models are in most cases larger than differences in the individual models across experiments. Therefore, our results highlight that the representation of transport, radiative transfer and microphysics, as well as the treatment/characterization of aerosol size distributions, plays a crucial role in the emergence of the widespread simulated temperature response. Since these processes interact with one another, it is not trivial to extract the exact reason of the differences among models. Below, we discuss the role of transport, radiation scheme, and the usage of interactive aerosol models and aerosol size distribution.

We first analysed the evolution of the mass concentration of SO₄ resulting from the volcanic plume, which acts as a tracer for transport affected only by formation and removal processes, and simultaneously represents the LW forcing for particles of an effective radius smaller than approximately 0.5 µm. All models following the ISA-MIP HErSEA protocol show a transport bias, i.e. preferred transport into the NH, which is in disagreement with REMAP-GloSSAC-2023 observations. This is in agreement with the findings of Quaglia et al. (2023), who relate the transport bias on the one hand to the point-injection strategy (compare also in our manuscript Figs. 2 and A3) and on the other to the absence of the Cerro Hudson eruption of August 1991 in the ISA-MIP HErSEA protocol. Such biases could be at least partially mitigated in future MIPs through application of a regional injection strategy (Tilmes et al., 2023), consideration of volcanic ash emissions alongside sulfur (Niemeier et al., 2009, 2021), and consideration of co-occurring volcanic events.

We also identified marked differences across models in SO₄ concentration already in the first months after SO₂ injection concerning the vertical spread of the SO₄ plume, the rate of confinement in the tropical pipe vs. partial transport through the shallow branch of the BDC. We attribute these early differences to the combined effect of SO₂ conversion rates, transport, model vertical velocity and resolution, in broad agreement with model differences identified in Niemeier et al. (2020). In the subsequent months, the models continue to differ in terms of the way of transport. While some models (e.g. EMAC, ECHAM6-SALSA, SOCOL-AERv2) show preferred transport through the shallow branch, others (e.g. UM-UKCA, WACCM6-CARMA, MIROC-CHASER) mainly transport through the deep branch. Similar differences in transport have also been highlighted in multi-model studies with a focus on age of air, in CMIP6 (Abalos et al., 2021) and CCMI-2020 (Eichinger et al., 2019). Furthermore, recent research revealed effects of QBO representation on transport (e.g., Zhuo et al., 2024). Following the experimental protocol, all models matched the observed post-Pinatubo QBO transition, either by nudging or from re-initialising from a prior control run. Given this uniform QBO treatment and the small effect on temperature anomalies in the models analysed (see Sect. 2.3.2 and Fig. A2), we do not consider it a prime contributor to overall model spread.

In addition, some differences in vertical transport/sedimentation are apparent in the models. Here, MIROC-CHASER is an outlier in the model set, with limited sedimentation and a non-apparent subtropical transport barrier. This might be to some extent explained by the models' low horizontal resolution (T42) and model top (∼ 50 km), although other low-resolution models such as SOCOL-AERv2 manage to represent the transport quite well. To analyse if the identified differences in transport dominate the model temperature response, we compared the co-alignment of the peaks in SO₄ concentration and temperature anomaly. While at first order we find a good agreement between these peaks in the individual models, significant differences exist among the models in the magnitude of the temperature anomaly and, consequently, in their sensitivity to SO₄ concentrations. These might be explained by differences in radiative transfer schemes (see Table 1) and furthermore differences in the particle size distribution (Quaglia et al., 2023).

Since for most of the models the effective radius is smaller or around 0.5µm (shown in Quaglia et al., 2023) in the first year after the eruption, we assume that absorption is not substantially affected by particle size. We note for completeness that ULAQ-CCM and WACCM6-CARMA have substantially larger effective radii at lower altitudes. Temperature anomalies in these layers differ, however, strongly between these two models, and as information on aerosol size distribution is not available beyond the effective radius, further interpretation remains unfeasible. Another exception is ECHAM6-SALSA which shows only a weak warming signal in our analysis. Such a dampened temperature signal for ECHAM6-SALSA was also identified by Laakso et al. (2022) in comparison to simulations performed with ECHAM6-M7. The authors attribute this dampening to the efficient formation of new particles occurring at the expense of gaseous sulfuric acid condensing onto pre-existing particles. Interestingly, however, the effective radius in ECHAM6-SALSA was found to be comparable to observations and the other models (Quaglia et al., 2023). Furthermore, transport and sedimentation also depend on particle size and substantial differences in the deposition of volcanic sulfur from the 1815 eruption of Mt. Tambora have been identified in Marshall et al. (2018). Recently Tilmes et al. (2026) highlighted in a study focusing on SAI sensitivity that the aerosol burden and mass distribution can vary significantly in the same modelling framework depending on the aerosol microphysical module used (sectional or modal), with implications for radiative forcing and impacts. Taken together, this shows that the simulated size distribution is somewhat important for a correct representation of the observed temperature response. Therefore, a detailed analysis of the size distributions of the models and the ability of the effective radius to represent them is suggested for future work. For now we consider this beyond the scope of the presented study given that for the bulk of the models effective radii do not seem to be a main contributor to the differences identified. Another candidate to explain the spread in temperature response in the ISA-MIP models, is the use of interactive aerosol schemes. To test this hypothesis we compared the ISA-MIP multi-model distribution with those obtained from CCMI-2022 and CMIP6-AMIP, which use prescribed SADs. If indeed interactive aerosol treatment would be a key contributor to model spread, one would expect to obtain narrower multi-model distributions from CCMI-2022 and CMIP6-AMIP. Interestingly, this hypothesis has to be rejected, as no reduction in the multi-model spread is found. One important point to note is that all MIPs considered comprise different model families. And, while the individual multi-model temperature anomalies agree quite well with those of reanalyses, our study shows that agreement in the median occurs through error compensation when averaging across models. The analysis and the sensitivity tests for model selection presented here highlight that by far the most important factor driving both magnitude and spread of the multi-model distribution in temperature response to volcanic aerosol forcing is model choice.

The results of our multi-model comparison argue for caution given the large spread in aerosol burden and transport and resulting stratospheric temperature response – even in the case of a strict experimental protocol. Even more so, given our present-day understanding of the role of aerosol-induced stratospheric warming for dynamics and chemistry and subsequently coupling across atmospheric layers and potentially surface climate.

Given the importance of progress in understanding the root cause of differences in model response to stratospheric (volcanic) aerosol forcing, some future directions can be informed from the results presented and limitations identified here: (i) implementation of strict experimental protocols to ensure a robust baseline for comparison, while understanding that this alone does not allow for a full attribution of similarities and differences across participating models. (ii) supplementing the core experimental protocol with process-oriented MIPs to individually address uncertainties in core components, such as radiative transfer, transport and sedimentation. (iii) sensitivity tests to identify the impact of changes in model subroutines affecting core variables for aerosol-climate connections to ensure robustness within model families and generations. (iv) more than one ensemble member per model simulation to quantify intra-model spread. (v) strict MIP protocol ensuring delivery of data variables necessary to support all items listed above.

Appendix A

Figure A1 Temperature anomaly distribution of reanalysis data before (bright colours) and after MLR (dark colours) in the tropics between 70 and 20 hPa in the period from September to December 1991. All calculations use 1980–2000 as the base period.

Figure A2 Temperature anomalies in the Med-22km experiment for selected ISA-MIP models with QBO effects included or regressed out. For ULAQ-CCM, MIROC-CHASER and SOCOL-AERv2 temperature anomalies are derived by using the Control runs as reference. For ULAQ-CCM_QBO, MIROC-CHASER_QBO and SOCOL-AERv2_QBO the QBO signal was removed in experiment and Control, using the wind field at 30 and 50 hPa before anomalies were derived.

Figure A3 Temperature anomalies in experiment Med-22km for UM-UKCA runs with different injection strategies. UM-UKCA uses the injection strategy as defined in the ISA-MIP HErSEA protocols (SO₂ injection in a single grid cell close to the Mt. Pinatubo location (15° N, 120° E)) while in UM-UKCA^* SO₂ is injected at Mt. Pinatubo longitude in a latitude range between 0 and 15° N. REMAP-GloSSAC-2023 observations are provided for convenient reference in the top row.

Figure A4 Vertical time evolution of the temperature anomalies for all experiments and models in the tropics (20° S–20° N). In all panels: SO₄ concentrations are shown as contours (levels are 0.5, 1, 2, 5 µg m⁻³), dots indicate non-significant areas (p-value > 0.05) for reanalysis and models with the exception of ULAQ-CCM where only one ensemble member is available, and the thick black vertical line represents the Mt. Pinatubo eruption time.

Figure A5 Vertical zonal temperature anomaly for all experiments and models for the period September to December 1991. In all panels: SO₄ concentrations are shown as contours (levels are 0.5, 1, 2, 5 µg m⁻³), dots indicate non-significant areas (p-value >0.05) for reanalysis and models with the exception of ULAQ-CCM where only one ensemble member is available, and the vertical lines embed the tropical region (20° S–20° N).

Figure A6 Vertical profile of effective radius (a) and time evolution of stratospheric effective radius (µm) (b) in the tropics for Med-22km. Model results are compared to SAGE II following Quaglia et al. (2023).

Table A1 Available CMIP6-AMIP and CCMI-2022 models.

Table A2 RMSE for the time evolution of the observed (multi-reanalysis mean) and model-derived tropical mean temperature anomaly normalized to June 1991 for different pressure levels for individual models.