Validating a microphysical prognostic stratospheric aerosol implementation in E3SMv2 using observations after the Mount Pinatubo eruption

This paper describes the addition of a stratospheric prognostic aerosol (SPA) capability – developed with the goal of accurately simulating sulfate aerosol formation and evolution in the stratosphere – in the Department of Energy (DOE) Energy Exascale Earth System Model, version 2 (E3SMv2). The implementation includes changes to the four-mode Modal Aerosol Module microphysics in the stratosphere to allow for larger particle growth and more accurate stratospheric aerosol lifetime following the Pinatubo eruption. E3SMv2-SPA reasonably reproduces stratospheric aerosol lifetime, burden, aerosol optical depth, and top-of-atmosphere flux when compared to remote sensing observations. E3SMv2-SPA also has close agreement with the interactive chemistry–climate model CESM2-WACCM (Community Earth System Model version 2–Whole Atmosphere Community Climate Model) – which has a more complete chemical treatment – and the observationally constrained, prescribed volcanic aerosol treatment in E3SMv2. Global stratospheric aerosol size distributions identify the nucleation and growth of sulfate aerosol from volcanically injected SO₂ from both major and minor volcanic eruptions from 1991 to 1993. The modeled aerosol effective radius is consistently lower than satellite and in situ measurements (max differences of ∼ 30 %). Comparisons with in situ size distribution samples indicate that this simulated underestimation in both E3SMv2-SPA and CESM2-WACCM is due to overly small accumulation and coarse-mode aerosols 6–18 months post-eruption, with E3SMv2-SPA simulating ∼ 50 % of the coarse-mode geometric mean diameters of observations 11 months post-eruption. Effective radii from the models and observations are used to calculate offline scattering and absorption efficiencies to explore the implications of smaller simulated aerosol size for the Pinatubo climate impacts. Scattering efficiencies at wavelengths of peak solar irradiance (∼ 0.5 µm) are 10 %–80 % higher for daily samples in models relative to observations through 1993, suggesting higher diffuse radiation at the surface and a larger cooling effect in the models due to the smaller simulated aerosol; absorption efficiencies at the peak wavelengths of outgoing terrestrial radiation (∼ 10 µm) are 15 %–40 % lower for daily samples in models relative to observations, suggesting an underestimation in stratospheric heating in the models due to the smaller simulated aerosol. These potential biases are based on aerosol size alone and do not take into account differences in the aerosol number. The overall agreement of E3SMv2-SPA with observations and its similar performance to the well-validated CESM2-WACCM makes E3SMv2-SPA a viable alternative to simulating climate impacts from stratospheric sulfate aerosols.

1 Introduction

Explosive volcanic eruptions are a significant source of aerosol forcing given their propensity to inject gas and particulate matter into the stratosphere. These gases form long-lived aerosol that can be spread around the globe (Robock, 2000). The impacts of the stratospheric aerosol loading from these eruptions are wide ranging, as exemplified in observational and modeling findings following the Pinatubo eruption in 1991: ozone depletion (Hofmann et al., 1994; Solomon et al., 1993; Portmann et al., 1996), surface temperature decreases (Parker et al., 1996; Soden et al., 2002), lower-stratosphere temperature increases (Labitzke and McCormick, 1992), reduction in global precipitation (Gillett et al., 2004), lowering of global sea level (Church et al., 2005), increases in cirrus cloud cover (Liu and Penner, 2002; Wylie et al., 1994), and increased diffusivity of incoming radiation (Robock, 2000) with resultant impacts on the net primary productivity of plants (Gu et al., 2003; Proctor et al., 2018; Greenwald et al., 2006). The extent of these impacts is dependent upon characteristics of the eruption (e.g., magnitude; Marshall et al., 2019) and climate state (Zanchettin et al., 2022). The foundation of these physical impacts comprises the stratospheric aerosol microphysical properties and chemical interactions that occur after an explosive volcanic eruption. Being able to accurately simulate these aerosol microphysical and chemical reactions in Earth system models is important for improving the fidelity of simulations. It also enables climate attribution work that clarifies the role of the climatic state versus characteristics of the eruption on downstream impacts. As a first step toward this goal, this paper presents a validation of a prognostic volcanic aerosol implementation within the Department of Energy (DOE) Energy Exascale Earth System Model version 2 (E3SMv2) (Golaz et al., 2022; Wang et al., 2020) against observational data from the Pinatubo eruption. Furthermore, E3SMv2 is compared with version 2 of the Community Earth System Model (CESM2) (Danabasoglu et al., 2020) with the Whole Atmosphere Community Climate Model version 6 (WACCM6; Gettelman et al., 2019), which shares many similarities with E3SM in its aerosol microphysical parameterizations but has more advanced atmospheric chemistry. This is to help identify any performance issues associated with a simpler chemical treatment in E3SMv2 and to serve as further validation of our implementation.

The stratosphere contains a persistent background layer of aerosol between the tropopause and ∼ 10 hPa that consists primarily (i.e., ∼ 95 % by mass; SPARC, 2006) of sulfuric acid droplets or sulfate aerosol (Junge et al., 1961). These aerosol are sustained by the tropospheric influx of gas-phase sulfate precursors carbonyl sulfide (OCS) and sulfur dioxide (SO₂) – along with smaller amounts of sulfate aerosol – from anthropogenic emissions and small-to-moderate volcanic eruptions (Hamill et al., 1997; Vernier et al., 2011). Intermodel comparisons of background sulfur (S) mass attribute 0.32 ± 0.050 Tg to OCS, 0.012 ± 0.007 Tg to SO₂, and 0.156 ± 0.051 Tg to sulfate aerosol with large intermodal differences due to differing model chemistry, removal processes, and dynamically driven transport and stratosphere–troposphere exchange (Brodowsky et al., 2023). Large-magnitude volcanic eruptions of the last century have led to huge injections of sulfur-containing species – mainly in the form of SO₂ (Guo et al., 2004a) – into the stratosphere (∼ 2–20 Tg S; McCormick et al., 1995), causing the formation of more numerous larger particles compared to background conditions (Deshler, 2008). These volcanic particles, which have an e-folding decay time of ∼ 1 year (Barnes and Hofmann, 1997), scatter incoming solar radiation back to space and act to cool the Earth's surface by a few tenths of a degree Celsius (Robock and Mao, 1995). Recent interest in the intentional stratospheric injection of sulfate precursors to recreate this phenomenon as a means to counter anthropogenic climate change (i.e., geoengineering or stratospheric aerosol injection (SAI); Caldeira et al., 2013) has driven a variety of Earth system modeling studies that examine the various climate and stratospheric chemistry implications of SAI scenarios (Tilmes et al., 2009; Kravitz et al., 2012, 2015; Kleinschmitt et al., 2018; Visioni et al., 2018; Weisenstein et al., 2022; Visioni et al., 2022). Multi-model comparisons of past volcanic eruptions (Zanchettin et al., 2016; Marshall et al., 2018; Timmreck et al., 2018; Clyne et al., 2021; Zanchettin et al., 2022; Quaglia et al., 2023) further quantify the downstream climate impacts from natural stratospheric injection events, identifying where models differ and where improvements can be made based on the observed climate impacts.

When simulating large-magnitude explosive volcanic eruptions, some climate models use prescribed volcanic forcing datasets as a way to reduce computational demand and to avoid uncertainties in prognostic aerosol formation. These datasets can estimate forcing based on satellite data, ground-based retrievals, ice core records, and other volcanic evidence (Toohey et al., 2016). One such dataset is the Global Space-based Stratospheric Aerosol Climatology (GloSSAC), which prescribes aerosol properties from a compilation of satellite, airborne, and ground-based observations (Kovilakam et al., 2020; Thomason et al., 2018). While GloSSAC and other prescribed datasets provide an accessible approach for incorporating volcanic forcing and validating model performance, prescribed aerosol products have limitations. Limited dataset availability and/or spatial coverage necessitates data interpolation within the forcing dataset, which may not accurately represent the volcanic forcing in some regions (Kovilakam et al., 2020). Additionally, prescribed aerosols will not respond to the dynamic state in free-running coupled climate simulations as the database has been generated from the observed climatic conditions. This artificially constrains the volcanic forcing across ensembles of simulations and creates a disconnect between volcanic forcing and the actual atmospheric transport patterns, limiting the usability of these simulations for detection of an evolving impact and its attribution to the volcanic source. Another limitation is the lack of aerosol microphysical representation and evolution from a volcanic eruption, which ignores aerosol indirect effects on clouds and does not allow for model feedbacks on aerosol size and lifetime.

Prognostically modeling the formation and evolution of sulfate aerosol from sulfur dioxide (SO₂) injected into the stratosphere is an alternative, more complete approach for simulating volcanic eruptions, with a variety of methods for representing sulfate aerosol mass, size, and number. This approach can serve to recreate conditions where observations are lacking as well as helping elucidate microphysical processes that contribute to aerosol properties. Aerosol forcing is also more dynamic in prognostic simulations given that it is not tied to the spatial pattern of the prescribed forcing. This allows for the simulation of evolving aerosol forcings and feedbacks in fully coupled model simulations or ensemble sets. The simplest prognostic approach is to use a bulk aerosol treatment, which prescribes an aerosol size distribution to a predicted aerosol species mass. This was applied to the earliest multi-year simulations of Pinatubo run with the Hamburg climate model ECHAM4 (Timmreck et al., 1999a, b) and has been used recently to show the large impact that SO₂ injection height can have on volcanic mass burden and climate forcing (Gao et al., 2023). The most accurate approach to simulating aerosol properties is the sectional (or bin) approach, but this can be computationally limiting depending on the number of aerosol size bins used. English et al. (2013) coupled the sectional Community Aerosol and Radiation Model for Atmospheres (CARMA; Toon et al., 1988) with version 3 of the Whole Atmosphere Community Climate Model (WACCM3; Garcia et al., 2007), showing the value of a sectional model in simulating the large variation in aerosol mode size and width that occurs after Pinatubo and larger-magnitude eruptions. More recently, Tilmes et al. (2023) showed that coupling CARMA to WACCM6 better represents the largest aerosol sizes following Pinatubo than a parallel-running modal aerosol model. The modal aerosol approach represents aerosol size distributions by multiple evolving lognormal functions. While this method strikes a balance between bulk simplicity and sectional cost, a downside is its dependence on defined modal widths, which can greatly impact stratospheric aerosol removal rates following Pinatubo if not tuned to match the observed stratospheric conditions (Kokkola et al., 2018). A modal aerosol approach is used in a stratospheric prognostic aerosol treatment with interactive ozone chemistry (Mills et al., 2016), developed in the Community Earth System Model (CESM, version 1) (Hurrell et al., 2013) using version 4 of WACCM (Marsh et al., 2013). This model design has been used to identify the impacts of volcanic eruptions on stratospheric ozone (Ivy et al., 2017; Solomon et al., 2016), the importance of interactive chemistry on the representation of sulfate formation and distribution following Pinatubo (Mills et al., 2017), and effective strategies for geoengineering (Kravitz et al., 2017).

The accurate simulation of the stratospheric sulfate size distributions is important for simulating the climate impacts of volcanic eruptions: both the scattering of incoming shortwave energy (i.e., surface cooling) and the absorption of outgoing longwave energy (i.e., stratospheric warming) are related to aerosol size and number through Mie theory, and size and number in turn depend on the aerosol modeling approach (e.g., bulk, bin, or modal). The choice of aerosol representation plays a large role in accurately representing the range of and variation in volcanic size distributions (English et al., 2013; Tilmes et al., 2023) and can result in a more than 50 % difference in the effective radius (R_eff) (Laakso et al., 2022) depending on a bin versus modal configuration. The performance of modal aerosol models is also highly dependent on the choice of modal widths (Weisenstein et al., 2007; Kokkola et al., 2018). The aerosol nucleation parameterizations within climate models can also affect aerosol formation, with the binary homogeneous nucleation scheme used in many climate models (Vehkamäki et al., 2002) potentially overpredicting nucleation rates by up to 4 orders of magnitude (Yu et al., 2023), leading to underestimating aerosol R_eff and overestimating aerosol number. Furthermore, neglecting van der Waals attractive forces, which aid in the coagulation of smaller aerosol particles, may also contribute to a small bias in R_eff in some models (English et al., 2013; McGraw et al., 2024).

Model simulations of Pinatubo use a variety of injection parameters, empirically chosen to match observations. The choice of vertical injection heights (17 km – Stenchikov et al., 2021; 18–21 km – Sheng et al., 2015; 18–20 km – Mills et al., 2016; 21–23 km – Dhomse et al., 2020) spans the 18–25 km range estimated from observations (Guo et al., 2004b), and models that do not include short-lived volcanic ash scale their mass emissions to account for the rapid removal of SO₂ that condensed on ash following the Pinatubo eruption (Neely and Schmidt, 2016; Mills et al., 2016; Clyne et al., 2021). The choice of SO₂ injection height can have a larger impact on volcanic mass burden and climate forcing than injection mass or particle size (Gao et al., 2023). Lower-stratospheric injection heights for Pinatubo (19 km) across model simulations result in too rapid a northward transport of the plume compared to observations, leading to more rapid removal and shorter aerosol lifetimes (a coarser model vertical resolution can lead to a similar effect; Brodowsky et al., 2021) (Quaglia et al., 2023). This is related to tropical aerosol retention, which is correlated with larger R_eff and a longer global mean aerosol optical depth (AOD) e-folding time (Clyne et al., 2021). Volcanic ash has a short lifespan on the order of days, and many models neglect to include this feature in their Pinatubo simulations. However, inclusion of volcanic ash in model simulations leads to strong absorption of longwave and shortwave radiation shortly after the eruption, leading to local dynamics changes in the cloud vicinity (Niemeier and Timmreck, 2009) and lofting the volcanic plume, the latter of which increases the plume height, aerosol lifetime, AOD, and R_eff (Stenchikov et al., 2021; Abdelkader et al., 2023).

Rates of sulfate aerosol formation and growth are driven by the oxidation of SO₂ by hydroxyl radical (OH) to form H₂SO₄ (sulfuric acid gas). This reaction in the stratosphere is OH limited. Models with more advanced interactive chemistry represent this OH depletion, resulting in longer aerosol lifetimes than models with simpler chemical treatments and prescribed OH (Bekki, 1995; Mills et al., 2017; Clyne et al., 2021). Prescribed OH models tend to rapidly oxidize available SO₂, which has been attributed to peak sulfate burdens occurring 3 months earlier than for interactive chemistry models studying the 1815 Tambora eruption (Clyne et al., 2021). Water vapor is also another important reactant in the formation of sulfate aerosol, controlling the nucleation of sulfuric acid gas into sulfate aerosol and increasing the availability of OH through its interaction with excited oxygen (O(¹D), a product of ozonolysis) to form two OH molecules (H₂O +O(¹D) → 2OH) (Seinfeld and Pandis, 2006; LeGrande et al., 2016). When it is coinjected with SO₂ from volcanic eruptions, water can significantly increase OH concentrations and plume AOD (LeGrande et al., 2016; Stenchikov et al., 2021; Abdelkader et al., 2023).

Here we present a new stratospheric prognostic aerosol capability within E3SMv2 that modifies the microphysical treatment of stratospheric aerosol in the four-mode Modal Aerosol Module (MAM4; Liu et al., 2012, 2016) to enable simulation of the evolution of volcanic stratospheric aerosols and their properties. Similarly to Mills et al. (2016), we add a stratosphere-specific sulfate treatment to complement the preexisting MAM4 chemistry and physics (default MAM4 includes the oxidation of SO₂ to form sulfate aerosols, their further growth through condensation and coagulation into larger aerosol size modes, sedimentation of these aerosols, and removal via wet and dry deposition). This model parallels work by Ke et al. (2024) and Hu et al. (2024) on the five-mode Modal Aerosol Module (MAM5) that incorporates more complete sulfate chemistry and an additional volcanic sulfate mode in E3SMv2. The validation of our implementation presented here will support forthcoming detection and attribution studies of societally relevant climatic impacts from stratospheric aerosols in free-running coupled climate simulations with varying volcanic source characteristics. By enabling dynamical consistency between transport, aerosol distribution, microphysical properties, and eruption characteristics (e.g., impact magnitude, timing, and location), this modeling capability facilitates the development of multivariate and multi-step attribution studies sensitive to spatio-temporal evolution (Hegerl et al., 2010). As future studies with this model capability will be free-running, they also enable better differentiation of the role of the climatic state in the detected and attributed impact.

With this new aerosol capability we detail the particle evolution and examine how model representations of the aerosol size distributions are related to global and regional radiative impacts at the surface and in the stratosphere. We use both observations from the 1991 eruption of Pinatubo and CESM2-WACCM to validate the implementation and demonstrate that we can reasonably simulate the lifetime, burden, AOD, and top-of-atmosphere (TOA) flux perturbations of stratospheric sulfate without the computationally expensive, whole-atmosphere, comprehensive chemistry in WACCM. These comparisons are enabled by constraining the model atmospheres to reanalysis horizontal winds present at the time of eruption. Additionally, we highlight the importance of simulated aerosol size and number for shortwave and longwave radiative impacts through comparisons of aerosol R_eff and number distributions to in situ observations. Utilizing these effective radii in single-particle Mie scattering calculations, we explore how variations in modeled aerosol microphysics affect the volcanic aerosol impacts on diffuse and direct radiation at the surface as well as longwave absorption and heating rates in the stratosphere.

2 Models and simulations

In this work we modify E3SMv2 to include a stratospheric prognostic aerosol capability. We test stratospheric microphysical and chemical implementations in the context of the Pinatubo eruption and compare them to version 2 of the Community Earth System Model (CESM2) with the Whole Atmosphere Community Climate Model version 6 (WACCM6; Gettelman et al., 2019).

2.1 E3SMv2

Originally a branch of CESM1, E3SMv1 diverged with a focus on computational efficiency, scalability, vertical and horizontal resolution, and aerosol and cloud parameterizations, as well as more physically based biogeochemistry, river, and cryosphere models (Leung et al., 2020). Simulations in this study are run with E3SMv2, in which clouds are parameterized with an improved version of the Cloud Layers Unified by Binormals (CLUBB) scheme (Larson, 2017), cloud microphysics are simulated by a two-moment bulk microphysics parameterization (MG2; Gettelman and Morrison, 2015), and mixed-phase ice nucleation depends on aerosol type and concentration as well as temperature (Wang et al., 2014; Hoose et al., 2010). Aerosols are simulated with MAM4 (Liu et al., 2016). The explosive volcanic eruption treatment prescribes stratospheric volcanic light extinction from version 1 of the GloSSAC reanalysis dataset (Thomason et al., 2018), which is mainly derived from SAGE II (Stratospheric Aerosol and Gas Experiment II) satellite measurements (Sect. 3.1.1, 3.5.1) assuming a sulfate refractive index from Palmer and Williams (1975) and aerosol volume and surface area assumptions from Thomason et al. (2008). E3SM uses the Rapid Radiative Transfer Method for GCMs (RRTMG, where GCMs stands for general circulation models) (Iacono et al., 2008; Neale et al., 2012), a two-stream approximation for calculating multiple scattering in the atmosphere from gas- and condensed-phase (i.e., aerosol, liquid cloud droplets, cloud ice, and hydrometeors) optical properties. Atmospheric chemistry is represented with version 2 of the interactive stratospheric ozone (O₃) model (O3v2; Tang et al., 2021). This model uses linearized stratospheric chemistry (Linoz v2; Hsu and Prather, 2009), which calculates net O₃ production as a function of temperature, local O₃ concentration, and overhead column O₃.

2.1.1 Prognostic aerosol in E3SMv2-PA

In the default prognostic volcanic aerosol simulations (E3SMv2-PA), prescribed volcanic extinction is removed and the sulfate aerosol precursor, SO₂, is emitted in the stratosphere. The emitted SO₂ undergoes chemical reactions to form sulfate aerosol and condenses onto the surfaces of preexisting aerosols following the prognostic calculations of MAM4.

In E3SMv2-PA, default MAM4 size modes are employed. All aerosol species are represented by three size modes, shown here with 10th- and 90th-percentile global, annual average number distribution dry-diameter ranges from Liu et al. (2012): Aitken (0.015–0.053 µm), accumulation (0.058–0.27 µm), and coarse (0.80–3.65 µm). A fourth mode (i.e., primary carbon mode (0.039–0.13 µm)) represents freshly emitted black carbon and organic carbon from combustion, which then ages into the accumulation mode. Aerosol mass and number are used to define a modal geometric mean diameter (D_g). This calculated D_g, in conjunction with a fixed modal geometric standard deviation (σ_g), defines the modal number distribution. The modal distributions then evolve based on nucleation (aerosol formation), evaporation (aerosol size reduction), condensation and coagulation (aerosol size growth), and dry/wet deposition (aerosol removal). In E3SMv2-PA, the growth of accumulation- to coarse-mode aerosol is not included because the troposphere seldom has high enough aerosol mass concentrations to generate such large aerosols through condensation and coagulation. Thus, accumulation-mode D_g is allowed to increase until it reaches the upper modal threshold (D_g,high), whereupon the model increases the accumulation-mode number to maintain D_g,high until D_g begins decreasing.

2.1.2 Prognostic stratospheric aerosol in E3SMv2-SPA

Because MAM4 was designed to accurately represent tropospheric aerosol at their respective concentrations and emission fluxes, the sulfate formation from the massive stratospheric influx of SO₂ from the eruption of Pinatubo is not accurately represented in E3SMv2-PA. In addition to removing prescribed volcanic extinction – as in E3SMv2-PA – further modifications are made to MAM4 to create a prognostic stratospheric aerosol version of E3SMv2 (E3SMv2-SPA) that has improved stratospheric aerosol representation following the Pinatubo eruption. These improvements borrow heavily from changes made to version 1 of CESM (CESM1) with WACCM (Appendix B in Mills et al., 2016), which are present in the default MAM4 version in CESM2-WACCM6. The major modifications to MAM4 include (1) the transfer of aerosol mass and number from the accumulation to coarse mode to increase aerosol size and represent the rapid aerosol growth following the Pinatubo eruption and (2) adjustment of the coarse-mode and accumulation-mode σ_g and minimum and maximum geometric mean diameters to increase aerosol lifetime. We note that these changes make the E3SMv2-SPA modal widths and size cutoffs identical to those in CESM2-WACCM6. Additional steps were taken to tune E3SMv2-SPA following the change in the accumulation- and coarse-mode size properties in (2), which included tuning of dust and sea salt emissions to account for the increased coarse-mode lifetime in the model as well as recalculating modal optical properties in MAM4 to account for the changes in aerosol size limits and distribution widths.

Unlike typical tropospheric conditions, explosive volcanic eruptions into the stratosphere provide ample SO₂ mass to drive sulfate aerosol into the coarse mode. To represent this rapid growth and overall larger aerosol diameters in the stratosphere, an irreversible accumulation-mode number and mass transfer into the coarse mode is added to E3SMv2-SPA. The model calculates the mass and number of particles in the tail of the distribution above a specified size cutoff (D_g,cut) of 0.44 µm, transferring this overshooting number and volume into the coarse mode. The model prohibits transfer from the accumulation mode if D_g<0.166 µm and allows total transfer of the grid cell mass and number when D_g>0.47 µm. In CESM2-WACCM, this transfer is reversible in the stratosphere, with an aqueous sulfuric acid (H₂SO₄) equilibrium pressure that depends on temperature and relative humidity. We left this out of our implementation under the assumption that, at the low relative humidities and low temperatures characteristic of the stratosphere, the effects from this process would be minimal.

To improve stratospheric aerosol lifetime in E3SMv2-SPA, the default coarse-mode σ_g is reduced from 2.0 to 1.2. The default accumulation-mode σ_g is also reduced, from 1.8 to 1.6, which has a small effect on aerosol lifetime. Additional changes to the aerosol modes allow for overlap between the coarse and accumulation modes and include increasing the accumulation-mode D_g,high from 0.44 µm to 0.48 µm and decreasing the coarse-mode lower threshold (D_g,low) from 1.0 to 0.4 µm. Lastly, coarse-mode D_g,high is increased from 4.0 to 40 µm. See Table S1 in the Supplement for a summary of these changes.

Most of the above changes have little effect on the tropospheric aerosols, except for changes to the coarse mode, which leads to longer-lived coarse-mode aerosol due to a reduction in removal rates. To account for this, emissions of dust and sea salt are tuned such that a simulation with perpetual present-day forcing obtains a global average AOD (0.1617) and global average dust AOD (0.0281) comparable to present-day remote sensing observations of ∼ 0.17 (Lee and Chung, 2013) and 0.028–0.03 (Ridley et al., 2016), respectively. The modal aerosol optical parameterizations are also affected by changes to the prescribed mode σ_g, and the modal optical properties were recalculated with the above modifications using the Ghan and Zaveri (2007) offline code used to generate the original files for CESM and E3SM.

The tuning of coarse-mode aerosol does not appear to significantly affect global measures of the simulated tropospheric climate. Two fully coupled, 164-year historical simulations (1850–2014) were run with E3SMv2-SPA, initialized from years 50 and 100 of a 100-year pre-industrial spinup simulation and run with the Model for Prediction Across Scales – Ocean (MPAS-Ocean) (Golaz et al., 2022). These simulations show the total AOD (Fig. S1), 2 m surface temperatures (T_2m; Fig. S2), and global radiative balance (Fig. S3) that track the five-member E3SMv2 historical simulations with prescribed volcanic forcing from Phase 6 of the Coupled Model Intercomparison Project (CMIP6; Golaz et al., 2022). Differences in atmospheric modes of variability (e.g., El Niño–Southern Oscillation (ENSO); Trenberth, 1997) due to internal variability affect T_2m during the Pinatubo period (Figs. S2, S4), but interval variability would average out if a mean were taken over more ensemble members.

2.2 CESM2-WACCM6

The major CESM2 atmosphere model improvements from CESM1 are the inclusion of CLUBB, MG2 cloud microphysics, MAM4, and orographic wave drag parameterizations (Danabasoglu et al., 2020). These are the same in E3SMv2 (Golaz et al., 2022) with the exception of a stratospheric prognostic aerosol treatment in MAM4 in CESM2-WACCM6 (Mills et al., 2016, 2017) (see Sect. 2.1.2). WACCM6 includes updated atmospheric chemistry, aerosol microphysics, and gravity wave drag parameterizations from previous versions of WACCM. Atmospheric chemistry is treated comprehensively through the whole atmospheric column, representing key chemical species and reactions across the troposphere, stratosphere, mesosphere, and lower thermosphere (Gettelman et al., 2019). Within the stratosphere and mesosphere, WACCM6 explicitly calculates the net production and transport of 97 different chemical species described by nearly 300 reactions (Mills et al., 2017). This comprehensive chemical treatment is a key difference between CESM2-WACCM6 and E3SMv2 (including E3SMv2-PA and E3SMv2-SPA), as E3SMv2 prescribes from observationally derived climatologies of OH and other relevant chemical species in sulfur and ozone chemistry (Hsu and Prather, 2009).

2.3 Simulations

E3SMv2 simulations are run from 1990–1993 (1989 discarded for aerosol spinup). The model uses the horizontal and vertical resolution described in Golaz et al. (2022), with the dynamics run on a ∼ 110 km horizontal grid, the physics run on a coarsened ∼ 165 km horizontal grid, and both dynamics and physics using the same 72-layer vertical grid with a model top at approximately 0.1 hPa. Simulations have prescribed sea ice and sea surface temperatures (Taylor et al., 2000) and nudged column-resolved U and V winds to 6-hourly Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2), reanalysis data (Gelaro et al., 2017). CESM2-WACCM6 simulations are run over the same time period on a 0.95° × 1.25° grid over 88 pressure levels (model top at ∼ 4.5 × 10⁻⁶ hPa). In addition to nudging model U and V winds to 6-hourly MERRA-2 reanalysis data, the CESM2-WACCM6 simulations also nudge column-resolved model temperature to the reanalysis to constrain temperature-dependent stratospheric chemical reactions. Note that both E3SM and CESM2-WACCM6 have a variety of stratospheric dynamics biases (e.g., Gettelman et al., 2019) that are avoided here through atmospheric nudging. An upcoming publication on E3SM stratospheric processes details a variety of biases in E3SM that may impact free-running volcanic eruption modeling, including a weak-amplitude tropical quasi-biennial oscillation which oscillates too frequently and a weak Brewer–Dobson circulation (Christiane Jablonowski, personal communication, 2024).

For simulations that prognostically simulate volcanic aerosol formation (E3SMv2-PA, E3SMv2-SPA, CESM2-WACCM6), the SO₂ emissions for explosive volcanic eruptions are from VolcanEESMv3.11, a modified version of Neely and Schmidt (2016). The VolcanEE3SMv3.11 dataset contains estimates of SO₂ from volcanic eruptions on a 1.9 × 2.5° latitude-by-longitude grid, with 1 km altitude spacing from the surface to 30 km. In our period of interest (1991–1993), this includes the Pinatubo, Hudson, Spurr, and Lascar eruptions. SO₂ emissions are provided in molecules cm⁻³ s⁻¹, and all eruptions occur over a 6 h period. The modifications to Neely and Schmidt (2016) are described in Mills et al. (2016) and include a reduction in SO₂ emissions for eruptions of over 15 Tg of SO₂ by a factor of 0.55 to compensate for missing ice and ash removal processes. In the case of Pinatubo, while 18–19 Tg of SO₂ erupted in the atmosphere, only ∼ 10 Tg remained in the stratosphere 7–9 d after the eruption (Guo et al., 2004b). This rapid reduction in SO₂ corresponds to >99 % removal of volcanic ash mass (Guo et al., 2004a). Therefore, 10 Tg of SO₂ is emitted in this dataset for further chemical and microphysical evolution (Mills et al., 2016). The emission takes place between 18–20 km, at a single latitude–longitude grid cell (i.e., no spreading). For all simulations, the VolcanEESMv3.11 file was merged with the monthly CMIP6 SO₂ emissions for non-explosive volcanic sources and then remapped to 1 × 1°.

A prescribed volcanic forcing simulation (E3SMv2-presc) is run in addition to the prognostic volcanic aerosol simulations. This simulation uses the default prescribed forcing dataset in E3SMv2 (GloSSAC V1) and allows for an additional validation of prognostic aerosol model performance where observational data are lacking.

Table 1 provides a summary of some of the key model characteristics for the different sensitivity studies used in this work. These studies include E3SMv2 with prescribed volcanic forcing and no emission of volcanic SO₂ (E3SMv-presc), E3SMv2 with the default MAM4 prognostic aerosol treatment (i.e., no stratospheric aerosol modifications) and emission of volcanic SO₂ (E3SMv2-PA), E3SMv2 with the prognostic stratospheric aerosol modifications and emission of volcanic SO₂ (E3SMv2-SPA), and CESM2-WACCM6 with emission of volcanic SO₂ (hereon referred to as CESM2-WACCM).

Table 1 Model details for the simulations used within this study. All simulations are run for 5 years (1989–1993) with 1989 discarded for aerosol spinup. All E3SMv2 simulations are run with U+V winds nudged to MERRA-2 reanalysis data; CESM2-WACCM has U+V winds and temperature nudged to MERRA-2 reanalysis.

2.4 Effective radius, size distributions, and Mie scattering calculations

Aerosol size distributions in the model provide information about how MAM4 represents volcanic aerosol evolution and can also help explain changes in radiation balance in the Earth system. Here we calculate the effective radius (R_eff) and use this in single-particle Mie scattering to understand how changes in size affect diffuse–direct radiation at the surface and absorption of longwave radiation in the stratosphere. We also plot stratospheric size distributions for our model simulations to visualize the evolution of volcanic injection of SO₂ into sulfate. See Appendix A and B for more details on these calculations.

3 Observational datasets

In addition to comparison to CESM-WACCM, we also employ observational datasets of sulfate burden, AOD, TOA flux, atmospheric temperatures, and microphysical properties (R_eff and size distributions) to substantiate the performance of the prognostic aerosol capability implemented within E3SMv2. Details of the observational datasets are presented below.

3.1 Sulfate burden

3.1.1 HIRS

The High Resolution Infrared Radiation Sounder (HIRS; Baran and Foot, 1994) is an infrared radiometer that has measured surface reflectance at 19 different infrared channels (3.7–15 µm) and one solar channel (0.69 µm) from a variety of polar-orbiting NOAA platforms since 1978 (Borbas and Menzel, 2021). This study uses HIRS-derived aerosol mass loading over the period of May 1991 through October 1993 from Baran and Foot (1994), who used the difference between 8.3 µm (aerosol sensitive) and 12.5 µm (aerosol insensitive) channels to isolate the transmission through the volcanic plume. To infer aerosol mass loading, they assumed an average stratospheric sulfate aerosol composition of 75 % H₂SO₄ + 25 % H₂O by mass, particle size, and concentration from dustsonde measurements in July 1991 (Deshler et al., 1992) and a single-scattering albedo calculated from Mie theory by integrating over scattering and extinction coefficients from an assumed lognormal distribution of radius 0.35 µm and standard deviation 1.6. These data cover 80° N–80° S at 5° resolution with a systematic error of ∼ 10 % (±1.4 Tg aerosol) due to assumptions in processing and uncertainty in background concentration. Additional minimum and maximum aerosol composition bounds are included in this range, namely 59 %–77 % H₂SO₄.

3.1.2 SAGE-3λ

This work also uses stratospheric sulfate burden taken from the SAGE-3λ dataset compiled for CMIP6 (ftp://iacftp.ethz.ch/pub_read/luo/CMIP6/, last access: 12 January 2023) as reported in Quaglia et al. (2023). The SAGE-3λ dataset uses the Stratospheric Aerosol and Gas Experiment II (SAGE II; Sect. 3.5.1) wavelengths of 0.454, 0.525, and 1.024 µm, fitting the measured extinction at these wavelengths to a lognormal size distribution and estimating sulfate mass burden from the number density, mode radius, and width of the distribution (Revell et al., 2017)

3.2 Aerosol optical depth

AVHRR

The Advanced Very High Resolution Radiometer (AVHRR; Zhao et al., 2013; Heidinger et al., 2014) is a radiometer that measures surface reflectance in six spectral bands (0.63, 0.86, 1.6, 3.75, 11, and 12 µm) and has served as a meteorological imaging sensor on the NOAA polar-orbiting platforms since 1978. It has a 1.1 km spatial resolution and, during the Pinatubo eruption period of interest (1990–1994), had two to four global views per day (Heidinger et al., 2014). An offline radiative transfer model is used to determine lookup tables for AOD retrievals. The radiative transfer model assumes fine- and coarse-mode aerosol properties based on validation of AVHRR with the surface radiometer measurements from the Aerosol Robotic Network (AERONET; Zhao et al., 2002), and uncertainty in the AVHRR AOD is estimated at 11.3 % based on surface AERONET validation (Zhao, 2022). This work uses monthly, clear-sky AOD retrieved over oceans and regridded to 1° resolution. Retrievals are made from channel 1 (0.63 µm) and related to the radiatively equivalent AOD at 0.6 µm through a radiative transfer and surface–atmosphere model (Rao et al., 1989). Detection limits on AOD from the AVHRR range from a minimum of 0.01 to a maximum of 2 (Russell et al., 1996).

3.3 Top-of-atmosphere radiative flux

ERBS

The TOA global radiative flux at a 1° × 1° resolution is used from version 2 of the Diagnosing Earth's Energy Pathways in the Climate system (DEEP-C) merged data product drawing from the Earth Radiation Budget Satellite (ERBS) near-global (60° S–60° N) non-scanning instrument and other reanalysis and observational datasets (Allan et al., 2014). The ERBS instrument measures reflected shortwave radiation and total outgoing radiation, allowing for the separation of longwave radiative flux by subtraction (Minnis et al., 1993).

3.4 Atmospheric temperature profiles

3.4.1 MERRA-2

The Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2), is a reanalysis product that assimilates satellite, radiosonde, radar, ship, buoy, and aircraft observations into version 5.12.4 of the Goddard Earth Observing System (GEOS) atmospheric general circulation model (Rienecker et al., 2011; Gelaro et al., 2017). These data are produced on a 0.5° × 0.625° grid with 72 vertical levels from the surface to 0.01 hPa. MERRA-2 observations include atmospheric state (temperature, pressure, humidity), dynamics, precipitation, radiation, and ozone, with updated aerosol observations from AVHRR over the period 1979–2002 (Gelaro et al., 2017).

3.4.2 RICH-obs

Version 1.5.1 of the Radiosonde Innovation Composite Homogenization (RICH-obs) software package is a compiled global radiosonde dataset that is merged with the help of reanalysis climatologies and neighboring data temperature records dating back to 1958 (Haimberger et al., 2012, 2008). The data gaps in station data are identified by divergence from 40-year climatology in the European Centre for Medium-Range Weather Forecasts Reanalysis (ERA-40), and the interpolation of these gaps is estimated from time series of neighboring radiosonde measurements, making RICH-obs less affected by satellite observations in the reanalysis but potentially biased in remote regions due to interpolation errors. This dataset also consists of 32 ensemble members that span a variety of sensitivity parameters and thresholds for interpolating to nearby radiosonde time series (Haimberger et al., 2012).

3.5 Effective radius and size distributions

3.5.1 SAGE II

The SAGE II (Mauldin et al., 1985) instrument flew from October 1984 to August 2005 on the Earth Radiation Budget Satellite (ERBS), measuring light extinction through the atmospheric limb at seven channels from 0.385 to 1.02 µm. The global coverage is 80° S–80° N with 1 km vertical resolution. This work uses aerosol R_eff from SAGE II version 7 (Damadeo et al., 2013) over the tropics (20° S–20° N), limited to 21–27 km (50–20 hPa) due to sparse data at lower altitudes (Quaglia et al., 2023). The R_eff is derived from a combination of extinction inversion algorithms that make use of the extinction ratios between 0.525 and 1.02 µm and assume that the aerosols are spheres with a sulfate composition of 75 % H₂SO₄ + 25 % H₂O by mass (Damadeo et al., 2013).

3.5.2 UARS–SAGE II (Stenchikov et al., 1998)

This work uses column average R_eff derived from two Upper Atmosphere Research Satellite (UARS; Grainger et al., 1995; Lambert et al., 1997) instruments (the Improved Stratospheric and Mesospheric Sounder (ISAMS) and Cryogenic Limb Array Etalon Spectrometer (CLAES)) and SAGE II extinction data (Stenchikov et al., 1998). The UARS instruments ISAMS and CLAES are limb sounders, reporting aerosol extinction at −12.11 and 12.66 µm wavelengths, respectively, across altitude profiles above 100 hPa (Lambert et al., 1997). Vertical resolution is approximately 2.5 km, and horizontal resolution is about 4° (Taylor et al., 1994). These extinction values are then used to derive aerosol R_eff assuming volcanic aerosol size distribution parameters from balloon-borne measurements (Deshler et al., 1992, 1993) and a sulfate refractive index corresponding to 75 % H₂SO₄ aerosol mass composition (Grainger et al., 1995; Lambert et al., 1997). Another limb sounder, SAGE II (McCormick et al., 1995), reports extinction at the 1.02 µm wavelengths, which is then used to derive aerosol number density in the atmosphere. The height-resolved number density and R_eff are then used to calculate column average R_eff zonal means at 40° N and 7° S from 200–10 hPa, reported in Fig. 4 of Stenchikov et al. (1998).

3.5.3 WOPC

The balloon-borne University of Wyoming optical particle counter (WOPC; Deshler et al., 1993; Deshler, 1994, 2003) uses particle scattering of white light to measure particle counts across 8–12 channels that range in size from 0.15–2 µm (Kalnajs and Deshler, 2022). These particle counts are then fit to unimodal or bimodal size distributions such that they minimize the root mean square error in number concentration between the measured cumulative count and the integral of the size distribution (Deshler, 2003). This instrument was launched from 1989–2013 from a variety of locations across the globe, with the most continuous sampling in Laramie, Wyoming. Instrument uncertainties include ±10 % in concentration, ±10 % in aerosol radius, and ±40 % in the distribution moments (e.g., surface area, volume, and extinction) (Deshler, 2003; Deshler et al., 2019). In this work we reproduce WOPC aerosol size distributions and derived R_eff over Laramie, Wyoming.

4 Results

This section traces model performance across scales (global to microphysical), tying the latter to the former through single-particle optical properties. The global performance of E3SMv2-SPA stratospheric mass burden (Sect. 4.1) and AOD (Sect. 4.2) is determined through comparisons with remote sensing data. Climate impacts of Pinatubo are explored via TOA radiative flux (Sect. 4.3) and atmospheric temperature (Sect. 4.4) perturbations. Regional comparisons to remote and in situ observations of stratospheric R_eff (Sect. 4.5) identify a model bias toward smaller sizes, which is explored in more detail through analysis of aerosol size distributions (Sect. 4.6). Lastly, model data are compared to individual balloon-borne measurements of aerosol size distributions and R_eff, relating the effective single-particle scattering (Sect. 4.7) and absorption (Sect. 4.8) efficiencies to changes in direct–diffuse radiation at the surface (Sect. 4.7.1) and longwave heating in the stratosphere (Sect. 4.8.1).

4.1 Sulfur burden

Figure 1 shows the stratospheric mass burden of the sulfur component of sulfate aerosol in the different model sensitivity tests, the HIRS observational dataset (Baran and Foot, 1994), and the SAGE-3λ dataset (Revell et al., 2017). When compared to HIRS and SAGE-3λ, E3SMv2-SPA improves the modeled stratospheric aerosol burden over E3SMv2-PA, especially in the years following the Pinatubo eruption. The increased aerosol burden – and thus, aerosol lifetime – in the stratosphere is mainly due to our modifications to the coarse-mode σ_g in E3SMv2-SPA. While E3SMv2-PA reaches a similar peak in sulfate burden, the underestimated aerosol burden following Pinatubo in E3SMv2-PA is mainly caused by too wide an aerosol number distribution, causing fast sedimentation of the larger coarse-mode particles in the upper tail of the distribution. The E3SMv2-SPA tends to overestimate aerosol burden compared to HIRS and SAGE-3λ in the 6 months after Pinatubo but agrees well with the slow decay reported in observations during 1992. In the 4 months following Pinatubo, models agree best with HIRS over SAGE-3λ, likely due to saturation issues identified in SAGE II limb-occultation data (Russell et al., 1996; Sukhodolov et al., 2018; Quaglia et al., 2023). From 1992 onward, stratospheric mass burden in E3SMv2-SPA agrees the best with SAGE-3λ, which reports higher burdens in 1993 than HIRS. E3SMv2-SPA and WACCM are similar in atmosphere and aerosol treatments but have very different atmospheric chemistry, which seems to impact lifetime.

In CESM2-WACCM, the interactive hydroxyl radical (OH) treatment causes OH depletion in the vicinity of the plume as the oxidation of SO₂ to form sulfate aerosol depletes available OH and therefore limits the reaction rate (Mills et al., 2017). This is in contrast to E3SMv2, which assumes an OH climatology unaffected by the oxidation of SO₂. The result is faster depletion of SO₂ and higher initial sulfate concentrations in the stratosphere in E3SMv2 (Fig. S5). The difference in OH treatment can be seen in Fig. 1, marked by the faster increase in the stratospheric sulfate burden in E3SMv2-PA and E3SMv2-SPA. Another significant difference between the two chemical treatments is the presence of carbonyl sulfide (OCS) in CESM2-WACCM, which is a largely inert tropospheric chemical species that is oxidized and photolyzed when it enters the stratosphere, forming sulfate. This is a major contributor to non-volcanic stratospheric sulfate and will lead to higher pre-Pinatubo sulfate concentrations in the stratosphere in CESM2-WACCM than in E3SMv2 (Mills et al., 2016). In Fig. 1, the effect of OCS is shown in larger pre-Pinatubo stratospheric burdens in CESM2-WACCM.

Figure 1 shows that E3SMv2-SPA performs reasonably well when forming sulfate from SO₂ and simulating increased aerosol lifetimes. The following sections will address how well the model parameterizes the aerosol microphysical properties and their impacts on the global radiative balance.

Figure 1 Stratospheric sulfate burden – reported in teragrams (Tg) of the sulfur mass contribution – for model simulations, as well as HIRS and SAGE-3λ remote sensing observations. The model data are processed to match the HIRS and SAGE-3λ data coverage of 80° N–80° S above the model lapse rate tropopause height. The sulfur component is determined by scaling modeled sulfate mass by the ratio of sulfur and sulfate molecular weights (MWs) such that Tg S = Tg SO₄ × $\frac{\mathrm{32.066}\mathrm{g}{\mathrm{mol}}^{-\mathrm{1}}}{\text{MW sulfate}}$. In both E3SMv2 and CESM2-WACCM, sulfate density and MW are assumed to be ammonium bisulfate ((NH₄)HSO₄; density = 1.7 kg m⁻³; MW = 115.11 g mol⁻¹) (Liu et al., 2012, 2016; Mills et al., 2016). Gray shading around the HIRS data represents systematic error of ∼ 10 % (±1.4 Tg aerosol) and the minimum and maximum aerosol composition bounds (59 %–77 % H₂SO₄).

4.2 Aerosol optical depth

In Fig. 2, the stratospheric contribution is isolated by subtracting the monthly mean AOD from pre-Pinatubo years, relying on the assumption that non-volcanic background aerosol in the atmosphere is similar in the near term. Model data are normalized to 1990 monthly means, while AVHRR is normalized to the monthly means for the period June 1989 to May 1991, with the exception of missing data from October–December 1990. Normalized model data are masked to reflect the same temporal and spatial sampling as AVHRR data from 1991–1993 over the oceans only and between 60° N–60° S. Here, AOD from the models is reported at 0.55 µm, while the AVHRR AOD is 0.6 µm.

The AVHRR AOD peaks 2 months after Pinatubo, linearly decreasing except for periods of flattening in the months March–July 1992 and January–July 1993. The flattening in 1992 is attributed to the continual growth of aerosols due to coagulation, which increases sulfate aerosol scattering efficiency (Russell et al., 1996; Stenchikov et al., 1998) even as aerosol burden continues to decay (Fig. 1). This is supported by a simulated increasing accumulation-mode – and, to a lesser extent, coarse-mode – D_g (Sect. 4.6; Fig. 7) and a slight increase in global average stratospheric aerosol R_eff (Sect. 4.5; Fig. 6). The flattening in 1993 for both observations and models may be due in part to the influence of the smaller Chilean volcanic eruption Lascar (30 January 1993; 23.36° S), which has a discernable impact on modeled global accumulation-mode aerosol mean diameters >0.2 µm (Sect. 4.6; Fig. 7).

In the first year following Pinatubo, AVHRR and the volcanically parameterized models (i.e., E3SMv2-SPA, E3SMv2-presc, and CESM2-WACCM) follow a similar trend, decreasing and leveling off near the beginning of 1992. As in aerosol mass burden, the overly short stratospheric aerosol lifetime in E3SMv2 leads to a rapid decay in AOD. All models consistently underpredict AOD in the first year after the eruption but tend to overpredict AOD in the third year following the eruption. Underprediction may be due in part to a lack of volcanic ash and incorrect number/size representation in the models. Overprediction is likely a factor of aerosol removal assumptions. Of the models, CESM2-WACCM has the closest agreement with AVHRR in 1993, due to a faster decline in AOD than E3SMv2-SPA. This indicates that the comprehensive chemistry in CESM2-WACCM may better represent aerosol size distributions than E3SMv2-SPA (see Sect. 4.4).

A surprising finding is the close agreement between E3SMv2-SPA and E3SMv2-presc. These models utilize very different approaches: prognostic aerosol microphysics from emitted SO₂ in MAM4 (E3SMv2-SPA) and prescribed stratospheric aerosol extinction from a range of observations in GloSSAC (E3SMv2-presc). As noted in Quaglia et al. (2023), the dependence of GloSSAC on SAGE II data – which are saturated at AOD ∼ 0.15 – means that E3SMv2-presc values may be underestimated in the months shortly after the eruption, suggesting that the close agreement between these datasets may be partly coincidental. Regardless, their similarities indicate that the changes to aerosol microphysics in MAM4 reasonably recreate the reanalysis data product, especially from 1992 onward when instrument saturation is less of a concern.

Figure 3 shows a zonal average of the stratospheric AOD values described in Fig. 2. As in the 60° S–60° N average, AVHRR exhibits much higher maximum AOD values in the stratosphere (0.428) than all other models (E3SMv2-PA: 0.175; E3SMv2-SPA: 0.204; E3SMv2-presc: 0.152; CESM2-WACCM: 0.25). The stratospheric parameterized models (E3SMv2-SPA (Fig. 3b), CESM2-WACCM (Fig. 3d)) have the closest agreement with AVHRR in the tropical plume shortly after the eruption. E3SMv2-PA, E3SMv2-SPA, and CESM2-WACCM fail to simulate the tropical confinement of aerosol present in AVHRR and E3SMv2-presc, simulating peak AOD that occurs further north and a weaker Southern Hemisphere AOD signal up to 1993 compared to observations and prescribed forcing datasets. Quaglia et al. (2023) identified this behavior in a range of models simulating the Pinatubo eruption, attributing this feature to model resolution, vertical wind structure, and the vertical distribution of the model volcanic cloud shortly after the eruption. Specifically, during the atmospheric conditions under which Pinatubo occurred, aerosols at levels <∼20 km are transported north while aerosols at levels >∼20 km are more effectively confined to the tropics (McCormick and Veiga, 1992). While E3SMv2-PA, E3SMv2-SPA, and CESM2-WACCM have the same injection parameters, E3SMv2-SPA has slightly more southern transport than CESM2-WACCM, which may be related to higher concentrations of sulfate above 20 km than CESM2-WACCM (Fig. S6). We speculate that this disagreement may arise from interactive chemistry in CESM2-WACCM and its effect on sulfate nucleation, growth, and SO₂ lifetime given the similarities in aerosol transport (i.e., nudged atmospheric dynamics), vertical resolution aerosol microphysics, and injection between the two models. E3SMv2-presc (Fig. 3c) has a noticeable low bias in AOD over the Northern Hemisphere over 1992, and we believe this is related to our use of an older version of the GloSSAC dataset (Thomason et al., 2018) in which higher latitudes have a low bias in AOD attributed to linear interpolation of the SAGE II data (Kovilakam et al., 2020). This may contribute to the close agreement between E3SMv2-SPA and E3SMv2-presc pointed out earlier, likely indicating a higher E3SMv2-presc signal with an updated GloSSAC dataset. Lastly, when compared to other model simulations in Quaglia et al. (2023), the volcanic Pinatubo signals in E3SMv2-SPA and CESM2-WACCM without land and AVHRR-missing-data masks (Fig. S7) show qualitatively similar patterns and magnitudes to ECHAM6-SALSA and ECHAM6-HAM with a comparable injection treatment (18–20 km; 7 Tg SO₂).

Figure 2 Stratospheric aerosol optical depth (AOD) over the ocean and across latitudes 60° S–60° N from the model simulations and AVHRR. Models calculate AOD at 0.55 µm, and AVHRR channel-1 AOD (0.63 µm) is processed to 0.6 µm. Both datasets are normalized to volcanically quiescent periods as described in the text. The Pinatubo eruption is marked with the dotted gray line at June 1991. The gray shading indicates ±11.3 % uncertainty in AVHRR AOD.

Figure 3 Stratospheric aerosol optical depth (AOD) zonally averaged over the ocean and across latitudes 60° S–60° N from the model simulations and AVHRR. Models calculate AOD at 0.55 µm, and AVHRR channel-1 AOD (0.63 µm) is processed to 0.6 µm. Both datasets are normalized to volcanically quiescent periods as described in the text.

4.3 TOA radiation flux

Figure 4 compares the global TOA radiative flux from model simulations to the all-sky ERBS observations over 1991–1993, subtracting out corresponding monthly means from the pre-Pinatubo year 1990 (note that this is a different non-volcanic period than used in previous publications (2001–2005 – Allan et al., 2014, and Liu et al., 2015; 1999 – Mills et al., 2017), which will result in differing magnitudes for ERBS over 1991–1993). Model TOA flux is shown for all-sky (solid lines), clear-sky (faint dashed lines), and aerosol-impact-only (faint dotted lines) conditions. The radiative flux is reported as absorbed shortwave radiation (ASR, positive downward flux; Fig. 4a), outgoing longwave radiation (OLR, positive upward flux; Fig. 4b), and net radiative flux (NET, positive downward flux; Fig. 4c). In Fig. 4a, ASR shows the clearest model separation 3–4 months after Pinatubo, corresponding to peak AOD (Fig. 2). There is close agreement between E3SMv2-SPA, E3SMv2-presc, and CESM2-WACCM during the year 1992, which corresponds to the largest sulfate particles during the Pinatubo plume evolution (see Sect. 4.5 and 4.6). The all-sky signal exhibits noise due to differences in atmospheric conditions (i.e., cloud cover, tropospheric aerosol) and surface albedo between the period of interest and our control year (1990). There is a clear seasonal increase in ASR in the 1991/92 and 1992/93 Northern Hemisphere winters relative to Northern Hemisphere summer. When clear-sky conditions (no influence from clouds) are compared to all-sky conditions in the models, the seasonality disappears, implying that the seasonality is cloud-related and cloud albedo was greater in the Northern Hemisphere winter of 1990 than the Northern Hemisphere winters of 1991/92 and 1992/93. Even with noise introduced by non-Pinatubo factors, there is a distinct all-sky ASR signal in E3SMv2-SPA, CESM2-WACCM, and E3SMv2-presc that is improved compared to ERBS.

The all-sky OLR (Fig. 4b), which is affected by both aerosol absorption of infrared emissions from the Earth's surface and the cooling of the troposphere and surface by the scattering of solar radiation, has a weaker response across these models than ASR. This is due in part to a less efficient absorption of outgoing longwave radiation than scattering of incoming solar radiation, leading to a lower sensitivity of OLR to aerosol growth and evolution (see Sect. 4.8). The largest spread in model simulations occurs during 1992 when aerosols are at their largest (i.e., highest absorption efficiency of longwave radiation; Sect. 4.8) and the highest reduction in surface temperatures was observed (Parker et al., 1996). All-sky E3SMv2-SPA has the greatest reduction in OLR from April 1992 to the end of 1993 and overestimates the longwave flux reduction compared to ERBS. This corresponds to E3SMv2-SPA overestimation of global AOD values compared to AVHRR over this period (Fig. 2). During this same period, CESM2-WACCM has slightly better agreement with ERBS, which may be related to the temperature nudging in this simulation which will modulate CESM2-WACCM surface temperature reduction and stratospheric temperature. When clear-sky OLR fluxes are compared, there is a weaker reduction in OLR for E3SMv2-PA, E3SMv2-SPA, and CESM2-WACCM and nearly no change in E3SMv2-presc during 1992. Due to the lack of stratospheric aerosol in E3SMv2-presc, this appears to be evidence of volcanic influence on high-altitude clouds which act to reduce OLR further, supporting conclusions from Liu and Penner (2002) and Wylie et al. (1994). Lastly, the aerosol-only model simulations remove the 1991/92 and 1992/93 wintertime peaks in the OLR signal, indicating similar or smaller OLR in 1990 compared to our period of interest due to cooler surface conditions.

The improvements in all-sky NET (Fig. 4c; solid lines) with volcanic parameterizations are less apparent across the models than in ASR (Fig. 4a) but do show improvement during the first 6 months after the eruption and during 1992. Differences in cloud cover and surface conditions between our period of interest and 1990 introduce substantial noise into this comparison, but the removal of clouds (clear-sky) and the isolation of aerosol TOA forcing (aerosol only) show a clear separation of volcanic parameterizing models and E3SMv2-PA.

Figure 4 Top-of-atmosphere radiative flux from model simulations and ERBS observations (Allan et al., 2014; Liu et al., 2015). The panels describe (a) absorbed solar radiation (ASR; positive downward flux), (b) outgoing longwave radiation (OLR; positive upward flux), and (c) net radiative flux (NET = ASR − OLR; positive downward flux). ERBS TOA flux is under all-sky conditions, while model TOA flux is shown under all-sky (solid line), clear-sky (faint dashed line), and aerosol-only (faint dotted line) conditions. Monthly mean data are normalized to the pre-Pinatubo conditions by subtracting respective monthly means from the year 1990, and both datasets are averaged over the entire globe.

4.4 Atmospheric temperature profiles

The radiation interactions described in Sect. 4.3 will lead to changes in atmospheric temperature, namely a warming of the stratosphere due to aerosol absorption of outgoing longwave radiation and a cooling of the surface due to reflection and scattering of incoming solar radiation by the aerosol plume. Figure 5 shows the 1992 annual mean atmospheric temperature anomalies (subtracting the 1990 annual mean) in the models (Fig. 5a–d), MERRA-2 reanalysis data (Fig. 5e), and RICH-obs radiosonde product (Fig. 5f). The year 1992 was chosen given the highest model spread in TOA flux (Fig. 4), peak modeled reduction in ASR (Fig. 4a) and reduction in OLR (Fig. 4b), and peak surface cooling (Parker et al., 1996) over this period. Models and observations share similar anomaly spatial patterns, with the exception of RICH-obs in the 60–90° S upper troposphere and near the tropical tropopause. Differences in RICH-obs may be related to temperature interpolation errors introduced in these remote regions due to fewer radiosonde datasets (Haimberger et al., 2012; Free and Lanzante, 2009). There is greater stratospheric warming in E3SMv2-SPA (Fig. 5b), E3SMv2-presc (Fig. 5c), and CESM2-WACCM (Fig. 5d) compared to E3SMv2-PA (Fig. 5a). Furthermore, there is an improvement in midlatitude warming at higher altitudes (i.e., 50 hPa and above) over E3SMv2-PA when comparing Fig. 5a–d to observations (Fig. 5f), reflecting the higher plume heights in these models (Fig. S6). CESM2-WACCM and MERRA-2 have very similar temperature magnitudes and distributions, which is due to temperature nudging of CESM2-WACCM to the latter reanalysis dataset. There is not as obvious a surface cooling difference between E3SMv2-PA and other models and observations. All datasets show a large cooling signal in the northern troposphere that roughly corresponds with early-1992 max AOD between 30 and 50° N (Fig. 3), but this cooling signal could be influenced by internal variability in the normalization year of 1990 (Sect. 4.3).

Table 2 shows 50 and 850 hPa pressure level averages from Fig. 5. These comparisons represent stratospheric and near-surface changes in temperature, with the 850 hPa level chosen to accommodate the lowest pressure level in the RICH-obs data. These latitude-weighted averages range from 65° S–65° N to avoid missing data in the upper atmosphere and surface RICH-obs data (Fig. 5f). This comparison shows stratospheric warming that is overestimated in E3SMv2-SPA (1.57 K) and underestimated in CESM2-WACCM (0.9 K) compared to MERRA-2 reanalysis (0.89 K) and previously reported estimates of ∼ 1 K (Ramachandran et al., 2000). RICH-obs struggles to represent lower-stratospheric warming due to either the aforementioned sparsity of data or their low horizontal resolution (5°) compared to models (1°) and MERRA-2 (0.5°) or both. E3SMv2-presc shows a stratospheric warming more than 3 times that of MERRA-2, which is likely due to a known error converting CLAES infrared extinction to the SAGE II- and GloSSAC V1-reported 1020 nm extinction coefficient, resulting in an exaggeration of peak aerosol extinction (Kovilakam et al., 2020). The 850 hPa cooling in CESM2-WACCM (−0.33 K) agrees best with MERRA-2 (−0.36 K) and RICH-obs ($-\mathrm{0.29}\pm \mathrm{0.007}$ K) anomalies, due in part to nudging of CESM2-WACCM temperatures to MERRA-2. There is small improvement in E3SMv2-SPA (−0.23 K) and E3SMv2-presc (−0.26 K) compared to E3SMv2-PA (−0.22 K), but it is unclear how much internal variability is influencing these values.

This comparison gives an all-sky snapshot of surface and stratospheric 1992 temperature anomalies due to Pinatubo. The 50 hPa height shows a clearer improvement in the simulated temperature anomaly in E3SMv2-SPA and CESM2-WACCM than the 850 hPa height due to the influence of interannual differences in internal variability (Sect. 4.3) and internal modes of variability (e.g., ENSO; Santer et al., 2014) in the troposphere. The model trends in stratospheric and near-surface temperature changes are consistent with changes in OLR and ASR (Fig. 4), respectively. Temperature trends also tend to agree better with observations and reanalysis with stratospheric volcanic parameterizations (E3SMv2-SPA, CESM2-WACCM) and prescribed volcanic aerosol (E3SMv2-presc). The next sections explore the microphysical representation within the models and how this influences lifetime, AOD, TOA flux, and temperature.

Figure 5 Annual mean change in atmospheric temperatures (K) for the year 1992, shown for (a–d) model simulations, (e) MERRA-2 reanalysis, and (f) the RICH-obs radiosonde data product. Anomalies are calculated by taking the difference between year 1992 and 1990 annual means. Model data are remapped from hybrid vertical coordinates to MERRA-2 pressure levels.

Table 2 Annual mean temperature anomalies at 50 and 850 hPa levels, shown from model simulations, MERRA-2 reanalysis, and the RICH-obs radiosonde data product. As in Fig. 5, anomalies are calculated as the difference between the year 1992 and 1990 annual means. Data are averaged over 65° S–65° N to avoid missing data in RICH-obs in the Antarctic. Included in RICH-obs is 1 standard deviation about the 32-member ensemble spread.

4.5 Stratospheric effective radius

R_eff has frequently been used to characterize stratospheric aerosol properties, with stratospheric R_eff of less than ∼ 2 µm leading to a net solar radiation scattering effect and surface cooling (Lacis et al., 1992; McGraw et al., 2024). Based on a range of in situ and remote sensing datasets, background stratospheric R_eff is estimated at 0.17–0.19 µm and following Pinatubo reaches average values of around 0.5 µm with observed values as large as 0.8 or 1.0 µm (Russell et al., 1996, and references therein). In the month following Pinatubo, there is little change in R_eff. This is due to a rapid increase in very small (i.e., Aitken-mode sulfate) and very large (i.e., ash) aerosol particles following the eruption. Enhanced coagulation and condensation, coupled with low sedimentation rates, lead to a steady increase in R_eff over the next 3–6 months. Aerosol growth continues until approximately mid-1992, when R_eff peaks, lagging peak values in other metrics such as mass burden and AOD. The smaller-magnitude eruptions of Hudson, Spurr, and Lascar also contribute to an increased R_eff over this period, with more of an impact in near-source regions.

Figure 6a shows global R_eff, in addition to R_eff from three different regional zones specific to different observational datasets: comparisons at 40° N and 7° S latitude bands and less than 100 hPa with UARS (Fig. 6b); 41° N, 105° W between 130–10 hPa with WOPC (Fig. 6c); and 20° S–20° N between 50–20 hPa with SAGE II (Fig. 6d). Observations tend to measure a minimum size that falls in the middle of the Aitken mode in the model. To account for this characteristic of the data, the model R_eff is an average of effective radii calculated with and without the Aitken mode, with the range between the two R_eff calculations represented by shading about the line. Maximum differences tend to occur before Pinatubo, shortly after Pinatubo, and with other volcanic eruptions (e.g., Hudson (8 August 1991; 45.9° S)). There is much larger spread between these two approaches in CESM2-WACCM in the Northern Hemisphere winter–spring, which may be due to enhanced stratosphere–tropopause exchange leading to higher concentrations of lower-stratospheric Aitken-mode particles (Sect. 4.6 and Fig. 7). In Fig. 6a the models reproduce the expected background R_eff of 0.17–0.19 µm, and the improvements to aerosol lifetime in E3SMv2-SPA and CESM2-WACCM can be seen in the slower decrease in R_eff compared to E3SMv2-PA. There is also a nearly identical pattern in E3SMv2-SPA and CESM2-WACCM data but with slightly higher R_eff in CESM2-WACCM.

All of these models underestimate R_eff compared to observations (Fig. 6b–d). In Fig. 6b, the finer temporal responses to the Pinatubo eruption in the UARS data are less apparent in the model, namely the large peak in R_eff at 7° S associated with short-lived volcanic ash in the stratosphere and the delayed peak at 40° N and 7° S due to particle aggregation (Stenchikov et al., 1998). The models neglect volcanic ash contributions, explaining the more gradual particle growth at 7° S. While the models do not show the same sensitivity to increases in R_eff at these latitude bands – possibly due to the higher vertical and horizontal spatial scale in the number concentrations from SAGE II observations – E3SMv2-SPA does show a flattening of the 7° S R_eff akin to the UARS estimate. This corresponds to the eruption of Hudson in Chile and the resulting high influx of smaller Aitken-mode particles into the southern stratosphere, which drives down R_eff. The sensitivity to this eruption in E3SMv2-SPA may be due to higher Aitken-mode production in this model than in CESM2-WACCM (see Sect. 4.7).

The tropical regions tend to have a larger R_eff than the midlatitudes in the simulations. This is true of UARS regions 6–12 months after the eruption and also the SAGE II (tropics; Fig. 6d) data. These larger R_eff values persist in the tropics a year or more after the eruption. While all models exhibit similar modal aerosol diameters (Fig. S8), the higher R_eff is correlated with higher number concentrations in all aerosol modes and a slower decrease in aerosol number (i.e., reduced sedimentation) (Fig. S9). The reduced removal is likely due to higher initial concentrations of SO₂ in the volcanic plume over the tropics (Fig. S10) contributing to more rapid local aerosol growth and a net positive aerosol production. Furthermore, the presence of the upwelling branch of the Brewer–Dobson circulation in this latitude band may help suspend larger aerosol species, slowing aerosol sedimentation rates and increasing their lifetime. The eruption of the Lascar volcano in northern Chile in around February 1993 also contributes to a bump in SO₂, as well as the Aitken- and accumulation-mode aerosol number at the 7° S band (Figs. S9, S10).

When comparing the models to WOPC (Fig. 6c) and SAGE II (Fig. 6d), E3SMv2-PA has the closest agreement to these datasets in its initial aerosol growth. This growth is more rapid than the other models and leads to a peak in R_eff that, while being closer to observed values, drops off precipitously. The R_eff values in E3SMv2-SPA and CESM2-WACCM have the best agreement with observational values and decay rate a year or more after Pinatubo. Differences across the models are due to the different microphysical assumptions, which can be explored by looking at aerosol size distributions.

Figure 6 Stratospheric aerosol effective radius (R_eff) averaged (a) globally above the model tropopause; (b) at 40° N and 7° S and <100 hPa to compare with UARS data (Stenchikov et al., 1998); (c) over Laramie, Wyoming (41° N, 105° W), at 130–10 hPa to match WOPC data (Quaglia et al., 2022); and (d) over 20° S–20° N and 50–20 hPa to match SAGE II observations (Quaglia et al., 2022). The shaded range in models represents R_eff with and without the Aitken mode, and the model line is an average of the two. Error bars in WOPC data assume a 40 % uncertainty and correlation coefficients of 0.5 between aerosol moments and levels (see Appendix A2 in Quaglia et al., 2023). The Pinatubo eruption is marked with the dotted gray line at June 1991.

4.6 Aerosol size distributions

While R_eff is a good representation of aerosol size in the context of optical properties, it is not always a good indicator for behavior of aerosol microphysical processes. For example, an increase in the accumulation-mode particle number and a decrease in the coarse-mode particle number could manifest as an unchanging R_eff. An examination of the aerosol size distribution can be more informative when understanding how aerosol chemistry and the MAM4 microphysics contribute to model performance.

Figure 7 shows the globally averaged stratospheric aerosol distributions from the three models used in this study and their evolution from 1991 to the end of 1993. The contour fill represents aerosol number ($\mathrm{d}N/\mathrm{dlog}D$; cm⁻³) with dashed and dotted lines indicating modal D_g (µm). In all three models, there is growth in the Aitken mode due to new sulfate particle formation following major and minor volcanic eruptions during this period: Pinatubo (June 1991; vertical dotted line), Hudson (August 1991), Spurr in the Aleutian Islands (27 June 1992; 61.3° N), and Lascar (February 1993). In E3SM (Fig. 7a–b), the prescribed concentrations of OH lead to a rapid oxidation of available SO₂ and higher concentrations of Aitken-mode aerosol compared to CESM2-WACCM (Fig. 7c), in turn resulting in higher aerosol number concentrations in E3SM. It also appears that CESM2-WACCM has higher tropospheric aerosol transport into the stratosphere based on a higher concentration of Aitken-mode aerosol across the period of interest. CESM-WACCM has seasonal peak concentrations occurring asynchronously with volcanic eruptions and corresponding to Northern Hemisphere winter (e.g., 1991 and 1992). The increased Aitken-mode number concentration is also seen at 40° N (Fig. S9). These peaks are attributed to a lower tropopause in Northern Hemisphere winter (i.e., enhanced troposphere–stratosphere exchange) and may also be due to the inclusion of OCS in CESM-WACCM.

The larger volcanic eruptions, Pinatubo and Hudson, inject enough SO₂ into the stratosphere for Aitken-mode aerosols to grow through condensation and coagulation into the accumulation mode, which then grow through coagulation into the coarse mode. The exception to this is in E3SMv2-PA, which lacks the ability to transfer sulfate mass into the coarse mode and so retains a roughly constant coarse-mode D_g derived from trace mass (10⁻²²–10⁻²¹ kg cm⁻³; global average) and number (10⁻⁵–10⁻⁴ cm⁻³; global average) concentrations of dust, sea salt, and sulfate aerosol advected from the troposphere (note that mass and number are related to aerosol size through Eqs. S1 and S2 in the Supplement). Aerosol growth through the modes can be seen in the increased number and an increasing trend in aerosol size. Dips in modal D_g correspond to sudden increases in aerosol number (e.g., Aitken-mode nucleation from freshly injected SO₂, transfer of a large aerosol number from Aitken to accumulation mode), while mass remains relatively unchanged. This leads to a division of mass across a larger number of aerosols and a subsequent decrease in D_g.

Overall, the aerosol modal diameters are similar across the three models, but slight differences exist. There are slightly larger accumulation-mode and coarse-mode D_g values in CESM2-WACCM (0.262 µm; 0.843 µm) compared to E3SMv2-SPA (0.259 µm; 0.749 µm). This could be because of the interactive chemistry in CESM2-WACCM, where lower nucleation rates lead to fewer Aitken-mode aerosols that initially grow faster through condensation as the longer-lived SO₂ in the stratosphere condenses on preexisting particles. Contrast this with E3SM, where high nucleation rates lead to more numerous smaller aerosols that consume the available SO₂ more quickly. The E3SMv2-PA model accumulation-mode D_g reaches a higher maximum D_g (0.285 µm) than E3SMv2-SPA or CESM2-WACCM (∼ 0.26 µm). This is due to the missing coarse-mode treatment in this model coupled with a larger D_g,high (Sect. 2.1.2). The result is an accumulation mode that grows to D_g,high following Pinatubo, whereby the number is increased to maintain this size instead of transferring the mass and number to a larger mode. This also explains the better initial agreement in R_eff between E3SMv2-PA and observations in Fig. 6b–d.

Figure 8 compares modeled aerosol size distributions to in situ measurements from WOPC, tracking the evolution of a slice of the plume from pre-Pinatubo (19 April 1991) to the end of 1993 (16 November 1993). Daily samples taken with the WOPC over Laramie, Wyoming, at a single level (18 km; roughly corresponding to the peak plume R_eff across the period (Fig. S11)) are used to validate daily data from E3SMv2-PA and E3SMv2-SPA. Modal D_g (vertical lines) and $\mathrm{d}N/\mathrm{dlog}D$ are denoted by dotted and dashed lines, as in Fig. 7. The effective diameter (D_eff) for models and WOPC is included to relate the changing distributions to aerosol–light interactions.

What follows is a breakdown of Fig. 8 into a rough timeline evolution in aerosol size distributions and D_eff.

Pre-Pinatubo conditions (Fig. 8a) are very similar between E3SMv2-PA and E3SMv2-SPA, and the D_eff values are nearly identical between model and observations. CESM2-WACCM has higher background number concentrations than E3SM.

A month after Pinatubo (Fig. 8b), all datasets exhibit rapid growth in the smaller-diameter aerosols, with some growth in the modeled accumulation mode but little sign of the coarse mode in E3SMv2-SPA. CESM2-WACCM has a massive increase in Aitken- and accumulation-mode number, possibly due to transport of SO₂ and Aitken-mode aerosol into the region from the tropics. This Aitken-mode peak is also seen at 17 and 19 km levels (Figs. S12–S13), resulting in a decrease in R_eff over this height range (Fig. S11), and may be related to lower-altitude sulfate aerosol in CESM2-WACCM (<20 km; Fig. S6) being more effectively transported into northern latitudes (McCormick and Veiga, 1992). The D_eff values are still similar across the different datasets, with the exception of CESM2-WACCM, which is smaller due to the Aitken-mode influence.

At 6 months after the eruption (Fig. 8c), a clear coarse-mode signal emerges in WOPC, E3SMv2-SPA, and CESM2-WACCM along with a sharp increase in the modeled accumulation-mode number. Models and WOPC are all comparable in their bimodal D_g and accumulation-mode number, while the WOPC coarse mode has a lower number than E3SMv2-SPA and CESM2-WACCM. E3SMv2-PA has a larger D_eff than other datasets (though it is still within the uncertainty in WOPC D_eff). This peak in D_eff – which is also noted in R_eff in Fig. 6 – is attributed to a wider accumulation mode and larger modal D_g in E3SMv2-PA. CESM2-WACCM has a larger modal number than all other datasets.

At 11 months after the eruption (Fig. 8d), the WOPC data continue to grow in the accumulation and coarse mode. The accumulation-mode number decreases while the coarse mode is increasing in E3SMv2-SPA, suggesting conversion of accumulation mass into the longer-lived coarse mode. The WOPC coarse-mode D_g (1.92 ± 0.19 µm) is more than twice that of E3SMv2-SPA (0.828 µm), while the WOPC accumulation-mode best fit (D_g=0.909 ± 0.09 µm, s_g=1.23) is close to the E3SMv2-SPA coarse mode (D_g=0.828 µm, s_g=1.2) and is nearly identical to the CESM2-WACCM coarse mode (D_g=0.918 µm, s_g=1.2). The overall larger aerosols in the WOPC distributions lead to the largest difference in D_eff between the model and WOPC.

For 1 year and longer after the eruption, the accumulation mode in E3SMv2-PA follows a steady and rapid decay that can be tracked in size and number (Fig. 8d–i). Over this period, E3SMv2-SPA and CESM2-WACCM generate very similar size distributions, with a slightly lower coarse-mode number and slightly larger coarse D_g in CESM2-WACCM. Overall, WOPC and E3SMv2-SPA/CESM2-WACCM have similar distributions and similar D_eff values, with E3SMv2-SPA showing a slightly smaller D_eff, but it is still within the uncertainty of WOPC data. A peak in the Aitken-mode aerosol in the models in Fig. 8h may be due to the Lascar eruption in February 1993.

The previous timeline shows an improved aerosol evolution in E3SMv2-SPA compared to WOPC and exhibits the similarities between CESM2-WACCM and E3SMv2-SPA. It also indicates that E3SMv2-SPA does not simulate large enough aerosol compared to the observations. Similar behavior in D_eff and relative size bias is noted at the 17, 19, and 20 km levels as well (Figs. S12–S13). Given that the model is nudged to meteorology and that the aerosols in the stratosphere tend to evolve relatively slowly over time, it is assumed that the resolution artifacts of comparing a 1° × 1° model grid cell to in situ measurements are minimal. We present a couple of possibilities for this size bias. One is that the aerosol density is too large due to the model assumption of ammonium bisulfate density of 1.7 kg m⁻³ instead of the sulfuric acid density of ∼ 1.6 kg m⁻³ (Seinfeld and Pandis, 2006). This may bias the aerosol to be small while also contributing to more rapid removal of coarse-mode aerosol. Another possibility is the lack of van der Waals forces in both E3SMv2 and CESM2-WACCM. This intermolecular attraction has been shown to aid in the coagulation of smaller particles, enhancing R_eff (English et al., 2013). The impact of this size bias on aerosol climate impact is now explored through effective single-particle scattering.

Figure 7 Globally averaged stratospheric lognormal aerosol size distributions from 1991–1993 for E3SMv2-PA, E3SMv2-SPA, and CESM2-WACCM. The dotted, dashed, and long-dash lines indicate geometric mean diameters of the Aitken, accumulation, and coarse aerosol modes, respectively. The color contour is lognormally scaled. The Pinatubo eruption is marked with the dotted gray line at June 1991.

Figure 8 Stratospheric aerosol size distributions from E3SMv2-PA (blue), E3SMv2-SPA (red), and WOPC (black) from 1991–1993. WOPC launches are chosen to correspond to spring, summer, and winter measurements for each year, and samples are taken from the 18 km measurements and matched to the nearest model height and grid cell over Laramie, Wyoming (41.3° N, 105° W). The dotted, dashed, and long-dash lines indicate geometric mean diameters (vertical) and $\mathrm{d}N/\mathrm{dlog}D$ size modes (curves) of the Aitken, accumulation, and coarse aerosol modes, respectively. Triangles denote effective diameters (D_eff) derived from the size distributions. Uncertainties in the WOPC diameter, number, and D_eff are denoted by gray bars at the peak of the distributions and on the markers. Note that the anomalously large maximum number in the CESM2-WACCM Aitken mode (376.5 cm⁻³) is out of the plotting bounds.

4.7 Single-particle scattering efficiency

The R_eff can be used to characterize the evolving aerosol size distribution, which conveniently can be used in the offline calculation of single-particle optical properties using Mie theory (see Appendix B). We choose to use offline Mie code for both the model and observations to allow for a direct comparison between optical properties derived from modeled and observed R_eff. Figure 9 shows the scattering efficiency (Q_s) of sulfate particles of size R_eff for the same time and height samples as Fig. 8. Additionally, this plot marks the approximate wavelength of peak solar blackbody irradiance (λ_solar=0.5 µm) to identify when aerosol scattering of solar radiation has the largest impact on radiative balance. Q_s is proportional to the effective size parameter ($x_{\mathrm{eff}}=\mathrm{2}\mathit{\pi }{R}_{\mathrm{eff}}/\mathit{\lambda }$) to the fourth power ($x_{\mathrm{eff}}^{\mathrm{4}}$) (i.e., λ⁻⁴) (Petty, 2006). This scattering feature can be seen in Fig. 9, where, as wavelengths increase, there is a rapid decrease in Q_s.

The peak in Q_s shifts from λ of ∼ 0.2 µm before and shortly after Pinatubo (Fig. 9a–b) to 0.5–0.7 µm 6 months after the eruption (Fig. 9c) with the onset of rapid particle growth. For the remainder of the time (Fig. 9d–i), E3SMv2-SPA and CESM2-WACCM R_eff result in a Q_s peak right around λ_solar, indicating a very efficient scattering of sunlight by these simulated aerosols. The WOPC reports larger R_eff than the models, which have a Q_s that is 10 %–80 % higher than observations at λ_solar during 1992 (Fig. 9d–f), becoming more similar during 1993 with the decay of stratospheric R_eff (Fig. 9g–i). E3SMv2-PA has a similar Q_s to E3SMv2-SPA/CESM2-WACCM until the end of 1992 (Fig. 9f), when there is a quick progressive drop in R_eff due to rapid stratospheric aerosol deposition.

These results indicate an optimal sulfate scattering R_eff of ∼ 0.3–0.4 µm, which is similar to estimates of the optimal sulfate aerosol size of 0.3 µm in geoengineering applications (Dykema et al., 2016). Given that WOPC predicts larger R_eff and smaller Q_s at λ_solar 6–18 months after Pinatubo, our results suggest that E3SMv2-SPA/CESM2-WACCM overestimates aerosol scattering and the global cooling effect of the Pinatubo aerosol at this level. This is supported by Fig. 4a, where both E3SMv2-SPA and CESM2-WACCM scatter more strongly than ERBS and E3SMv2-presc from April–August 1992 (though differences in clouds may also be contributing to this result). Temperature comparisons at 850 hPa over 1992 do not reflect a cooler surface in E3SMv2-SPA and CESM2-WACCM compared to observations (Table 2), but this comparison is not ideal for making a connection between perturbations in scattering and surface temperature due to the unknown role of internal variability on temperature. Figure 9 relates size distribution differences in models and WOPC to climate impact at a single level. For a better validation of the stratospheric plume and to see if the model disagreement changes when looking across multiple levels, the next step is to calculate Q_s from the stratospheric R_eff.

Figure 10 shows monthly Q_s generated from averaged stratospheric R_eff at a midlatitude site (solid line) (Fig. 6c) and the tropics (dotted line) (Fig. 6d). As in Fig. 9, a similar pattern is shown where smaller modeled R_eff values lead to higher modeled Q_s values at λ_solar compared to observed Q_s. However, the differences between modeled and observed Q_s at λ_solar are not as stark for the monthly stratospheric average values, with maximum differences on the order of 5 %–14 % as opposed to the 10 %–80 % seen in the daily 18 km level data. The larger R_eff in the tropics leads to a higher Q_s at λ_solar 1 month after Pinatubo (Fig. 10b) and lower Q_s at λ_solar 18 months after Pinatubo (Fig. 10f) compared to midlatitude measurements in the stratosphere. The E3SMv2-SPA and CESM-WACCM maximum Q_s values from both regions hover around λ_solar, with higher modeled Q_s at λ_solar in the tropics 2 years after the eruption (Fig. 10h) due to a consistently larger R_eff in this region.

Figure 9 Daily 18 km scattering efficiency (Q_s) using the effective size parameter ($x_{\mathrm{eff}}=\mathrm{2}\mathit{\pi }{R}_{\mathrm{eff}}/\mathit{\lambda }$) and sulfate refractive index at 0 % relative humidity (Hess et al., 1998). Effective radii (R_eff) from E3SMv2-PA, E3SMv2-SPA, and WOPC are calculated at 18 km and are marked by colored triangles at the top of the plot. The bottom x axis is wavelength (λ; 0.1–3 µm). The top axis is R_eff (0.1–1 µm) used in the calculation of Q_s via x_eff. The vertical turquoise line marks the solar blackbody wavelength of maximum irradiance (0.5 µm).

Figure 10 Monthly stratosphere scattering efficiency (Q_s) using the effective size parameter ($x_{\mathrm{eff}}=\mathrm{2}\mathit{\pi }{R}_{\mathrm{eff}}/\mathit{\lambda }$) and sulfate refractive index at 0 % relative humidity (Hess et al., 1998). Solid lines are scattering efficiencies over 41° N, 105° W, and dotted lines are scattering efficiencies over 20° S–20° N. Monthly R_eff values are from E3SMv2-PA, E3SMv2-SPA, CESM2-WACCM, WOPC, and SAGE II from Fig. 6c–d and are marked by colored triangles (41° N, 105° W; WOPC) and hollow diamonds (20° S–20° N; SAGE II) at the top of the plot. The bottom x axis is wavelength (λ; 0.1–3 µm). The top axis is R_eff (0.1–0.7 µm) used in the calculation of Q_s via x_eff. The vertical turquoise line marks the solar blackbody wavelength of maximum irradiance (0.5 µm).

4.7.1 Modeled direct–diffuse radiation

The differences in R_eff between the presented model and observational datasets ultimately affect diffuse and direct radiation at the surface by changing AOD and other aerosol optical properties. In scattering incoming shortwave radiation, a small fraction of incident light is scattered back to space by the volcanic aerosol, reducing the amount of energy incident to Earth. More substantially, the forward-scattered shortwave radiation increases the diffusivity of incident radiation. Crop yields are differentially impacted by the increase in diffuse radiation as certain canopy structures benefit from diffuse radiation penetrating to deeper leaves. Thus, while a decrease in direct radiation has an overall negative influence on crop yield, the increase in diffuse radiation allows some crop types to be less influenced (e.g., rice, soybean, and wheat) then others (e.g., maize) (Proctor et al., 2018). The Mauna Loa observatory documented a lessening of the direct beam and increase in the diffuse radiation after the eruption of Pinatubo (Robock, 2000), and globally there was a 21 % reduction in direct radiation and a 20 % increase in diffuse radiation at the surface (Proctor et al., 2018).

We diverge from the idealized Q_s calculations in Sect. 4.7 to analyze modeled radiation diagnostics in the atmosphere, calculated via the online model Mie code and a two-stream approximation for calculating multiple scattering in the atmosphere (Iacono et al., 2008; Neale et al., 2012). The latter calculates diffuse and direct radiation at the surface from incident radiation, AOD (i.e., the column-integrated β_s; Eq. B2), single-scattering albedo (SSA; Fig. S15), and aerosol asymmetry parameter (g; Fig. S16). In a general sense, the input aerosol optical properties determine what fraction of direct radiation makes it through the aerosol layer (i.e., AOD), how much of the incident radiation is absorbed in the aerosol layer (i.e., 1−SSA), and the degree to which the radiation is scattered by the aerosol layer (i.e., g).

Figure 11 shows the perturbation in diffuse and direct radiation at the surface due to the Pinatubo aerosol layer in the visible to ultraviolet wavelengths (with similar behavior in the near-infrared (Fig. S17)). The values are averaged over the whole Earth (Fig. 11a), the 41° N latitude band (Fig. 11b), and 20° S–20° N tropical region (Fig. 11c). There is an increase in diffuse radiation across all of the models and a mirrored decrease in direct radiation. These changes track the increasing R_eff (Fig. 6) and the enhanced Q_s(λ) (Fig. 10). The magnitude of the direct radiation loss tends to be larger than the increase in diffuse radiation at the surface because some of the scattered incoming solar radiation is lost to space. In all regions, E3SMv2-SPA and CESM2-WACCM have more diffuse and less direct radiation than E3SMv2-PA, the latter of which rapidly approaches zero as the aerosol are removed from the atmosphere. The peak direct and diffuse radiation influence occurs in the tropics (Fig. 11c) prior to the 41° N band (Fig. 11b) as aerosol transport to northern latitudes relies on stratospheric circulation. In the 41° N band (Fig. 11b) CESM2-WACCM has a higher proportion of diffuse radiation at the surface than all of the other models.

Diffuse radiation is not only attributed to aerosol size and material properties – as is Q_s in Fig. 10 – but also related to the aerosol number concentration. This is also true of AOD, which is a cloud-free, aerosol-specific input into multiple-scattering radiation calculations in the models. AOD, which accounts for both number and Q_s, gives a cleaner signal of the aerosol influence. Figure 12 shows AOD over the same spatial regions as Fig. 11. While peaks in AOD and diffuse radiation differ, the intermodal relationships are very similar between the two comparisons. As in Fig. 11b, Fig. 12b shows CESM2-WACCM with higher AOD over the 41° N band. This is attributed to both the larger aerosol R_eff (Fig. 6c) and higher accumulation- and coarse-mode number concentrations (Fig. S9a) in CESM2-WACCM. The higher number concentrations are correlated with more SO₂ transport into the midlatitudes in CESM2-WACCM (Fig. S10a), which is likely due to slower oxidation of SO₂ in CESM2-WACCM due to the aforementioned OH depletion in this model.

The impact of clouds on diffuse–direct radiation in the midlatitudes shows up in Fig. 11b, with peaks in normalized total cloud cover (right axis) corresponding to diffuse (direct) radiation peaks (troughs), especially in the late spring and early summer of 1992 and 1993. The approximate co-location of peaks in cloud cover (Fig. 11b) and AOD (Fig. 12b) in the midlatitudes in June 1992 and 1993 may indicate the presence of aerosol indirect effects. This is consistent with findings by Liu and Penner (2002) that show high rates of homogeneous ice nucleation persisting in the northern midlatitudes through July 1992.

Lastly, E3SMv2-presc has good agreement with E3SMv2-SPA and CESM2-WACCM when compared globally or in the tropics but reports ∼ 60 % of the modeled magnitudes for diffuse–direct radiation and ∼ 40 %–60 % of modeled AOD in the midlatitudes (Figs. 11 and 12, respectively). This may be due to a low bias in AOD in GloSSAC V1 (Thomason et al., 2018) used here where higher latitudes have a low bias in AOD attributed to linear interpolation of the SAGE II data (Kovilakam et al., 2020).

Figure 11 Pinatubo-induced changes in diffuse and direct solar radiation at the surface in the 0.2–0.7 µm band (left axis) and total cloud cover fraction (right axis; blue). Regions include global (a), midlatitude band (b; 41° N), and tropical band (c; 20° S–20° N) averages. The Pinatubo eruption is marked with the dotted gray line at June 1991, and the sample months for Fig. 10 are marked with the dash-dotted gray lines. The data are normalized to 1990 monthly means.

Figure 12 Pinatubo-induced changes in aerosol optical depth (AOD) at 0.55 µm. Regions include global (a), midlatitude band (b; 41° N), and tropical band (c; 20° S–20° N) averages. The Pinatubo eruption is marked with the dotted gray line at June 1991, and the sample months for Fig. 10 are marked with the dash-dotted gray lines. The data are normalized to 1990 monthly means.

4.8 Single-particle absorption efficiency

While the Pinatubo eruption resulted in a net climate cooling effect due to increased scattering, it also contributed to stratospheric warming through the absorption of outgoing longwave radiation within the aerosol layer (Kinne et al., 1992). Figure 13 shows the Q_a calculated offline from model and in situ R_eff, with a dark-red line indicating the wavelength of peak terrestrial blackbody irradiance (10 µm; λ_earth). While sulfate scatters strongly at visible wavelengths (λ_solar in Figs. 10 and 9), it is not an effective absorber in the visible band (Q_a=0 at λ_solar), and the determining factor in stratospheric heating is the magnitude of Q_a at λ_earth. The wavelength dependence of Q_a is weaker than Q_s (proportional to x_eff (i.e., λ⁻¹)) due to a weaker longwave absorption sensitivity to R_eff (Lacis, 2015), and Q_a is strongly tied to the imaginary part of the sulfate refractive index at 0 % relative humidity (n_Hess; Appendix B) (Hess et al., 1998). The proportionality to x_eff is denoted in Fig. 13 by the linear increase in Q_a(λ) with increasing R_eff; the dependence on the imaginary part of n_Hess is reflected in the unchanging pattern in absorption magnitudes for all datasets and times. Generally speaking, Q_s is characterized by changes in the wavelength and is more sensitive to aerosol size fluctuation (Figs. 10 and 9), while Q_a shows changes in the absorption efficiency and is less sensitive to aerosol size fluctuation.

In Fig. 13, the linear dependence of Q_a on R_eff means that an individual particle in the tropics has more absorption capacity than the midlatitudes across the whole period, and the observations will show more absorption than the models. This is explained by R_eff in these two regions shown at the top of each panel. For example, the smaller modeled R_eff at 41° N leads to a modeled Q_a that is 15 %–30 % lower than that of WOPC through 1992 (Fig. 13c–f).

Because the same refractive index is used across all models and observational datasets, the differences in Q_a are due only to differences in aerosol size. However, this is a simplifying assumption for the modeled Q_a values due to their explicit calculation of the soluble aerosol refractive index (n_s) in the stratosphere based on aerosol water and sulfate content (Eq. B4). E3SM simulates a stratosphere with unrealistically low-stratospheric water vapor (Keeble et al., 2021; Christiane Jablonowski, personal communication, 2024), which results in aerosols that are mostly sulfate (Fig. S18) with a refractive index similar to n_Hess (Appendix B). Higher water vapor in CESM2-WACCM means that the assumption of n_Hess in the calculation of Q_a does not account for the volume weighting of the hydrated aerosol refractive index (Appendix B). When interstitial stratospheric water and sulfate are used to calculate a volume-weighted n_s using n_Hess and n_wat, Q_abs at λ_earth is smaller than that reported with only n_Hess by ∼ 14 %–20 % over the 41° N band and by ∼ 13 %–18 % over 20° S–20° N (Table S3). This is due to a smaller-magnitude imaginary refractive index for water at this wavelength.

It is unclear whether injecting volcanic H₂O would help offset the stratospheric dry bias in E3SM. Abdelkader et al. (2023) showed that the injection of H₂O for a 20 km Pinatubo-like injection increases sulfate mass and stratospheric AOD by ∼ 5 %, but at the colder temperatures in the lower stratosphere almost all of this water vapor freezes and sediments out. Retention of water vapor in the stratosphere depends on injections at higher altitudes where temperatures are warmer. Given that E3SMv2-SPA already tends to produce higher sulfate burdens compared to other simulations (Fig. 1), it is possible that the enhanced sulfate aerosol mass attributed to the addition of water vapor in the plume would bias the aerosol sulfate content to be high compared to observations, even as it would improve upon the aerosol size. Furthermore, the lower injection height of our simulations (20 km) and a cold-point tropopause that is too cold in E3SMv2 (Christiane Jablonowski, personal communication, 2024) would aid in the rapid removal of most of the injected water, reducing the effect of the water injection on aerosol properties.

Figure 13 Monthly stratosphere absorption efficiency (Q_a) using the effective size parameter ($x_{\mathrm{eff}}=\mathrm{2}p{R}_{\mathrm{eff}}/\mathit{\lambda }$) and sulfate refractive index at 0 % relative humidity (Hess et al., 1998). Solid lines are scattering efficiencies over 41° N, 105° W, and dotted lines are scattering efficiencies over 20° S–20° N. Monthly R_eff values are from E3SMv2-PA, E3SMv2-SPA, CESM2-WACCM, WOPC, and SAGE II from Fig. 6c–d and are marked by colored triangles (41° N, 105° W; WOPC) and hollow diamonds (20° S–20° N; SAGE II) at the top of the plot. The bottom x axis is wavelength (λ; 3–40 µm). The top axis is R_eff (0.1–0.7 µm) used in the calculation of Q_a via x_eff. The dark-red line marks the terrestrial blackbody wavelength of maximum irradiance (10 µm).

4.8.1 Modeled longwave heating rates

The absorption of outgoing longwave radiation in the stratosphere results in local heating of the aerosol layer. Figure 14 shows the longwave heating rate (LWH; K d⁻¹) due only to Pinatubo aerosols and normalized to 1990 monthly means over the same regions as Figs. 11 and 12. As in Figs. 11 and 12, Fig. 14 is calculated from the online model radiation code. Higher initial global (Fig. 14a) and 41° N (Fig. 14b) LWH in CESM2-WACCM up to March 1992 mirrors model behavior for AOD (Fig. 12a–b). The similarity to AOD indicates a LWH dependence on aerosol absorption optical depth (AAOD; level-integrated β_a; Eq. B3) – and therefore aerosol number concentration and Q_a – in these two regions. As with diffuse–direct radiation (Sect. 4.7.1), characterizing the stratospheric longwave heating rate based solely on Q_a (i.e., aerosol size) is an oversimplification of the actual stratospheric absorption because number concentration is also acting to scale b_a (note that the E3SMv2-presc has no Pinatubo signal in LWH, which is related to how the prescribed volcanic forcing is treated in the model and is not indicative of a missing physical mechanism; as a result, E3SMv2-presc is not included in this comparison).

Figure 14 Pinatubo-induced changes in longwave heating rate (K d⁻¹) over the modeled longwave wavebands (waveband midpoints of ∼ 3.5–514 µm). Regions include global stratospheric mean (a), a midlatitude band stratospheric mean at WOPC pressure levels (b; 41° N, 130–10 hPa), and a tropical band stratospheric mean at SAGE II pressure levels (c; 20° S–20° N, 50–10 hPa). The Pinatubo eruption is marked with the dotted gray line at June 1991, and the sample months for Fig. 10 are marked with the dash-dotted gray lines. The data are normalized to 1990 monthly means, and contributions from non-aerosol effects (i.e., clouds and gases) have been removed.

5 Discussion

The preceding sections present encouraging results for E3SMv2-SPA and its ability to model explosive volcanic eruptions. These results begin at the global scale and move into the micro-scale, tying the latter to the former to understand the strengths and weaknesses of the modal aerosol modeling approach. While this stratospheric parameterization is new to E3SM, it is based on a modeling framework developed for CESM2-WACCM, which has an additional interactive chemistry element not present in E3SM (Mills et al., 2016). Along with a variety of observational datasets, E3SMv2-SPA is validated against CESM2-WACCM, which elucidates the effects of interactive chemistry on aerosol formation. Ongoing work on a five-mode version of MAM in E3SM (MAM5) will include interactive stratospheric sulfate chemistry and will have a separate stratospheric mode for volcanic sulfate to avoid impacts on tropospheric coarse-mode aerosol (Ke et al., 2024; Hu et al., 2024).

On the global scale, E3SMv2-SPA performs well compared to observational datasets and has similar behavior to CESM2-WACCM. The near globally averaged stratospheric sulfate burden in E3SMv2-SPA agrees well with SAGE-3λ and HIRS remote sensing datasets (Fig. 1), indicating that the model does a reasonable job of simulating aerosol lifetime and loading in the stratosphere. The higher mass burdens in E3SMv2-SPA compared to in CESM2-WACCM point to the OH consumption in CESM2-WACCM's interactive chemistry limiting aerosol formation from SO₂. For near-global AOD, the volcanically parameterized E3SMv2-SPA, CESM2-WACCM, and E3SMv2-presc all have similar performance and all exhibit a similar shape to the retrieved AVHRR AOD (Fig. 2). The models underestimate AOD compared to AVHRR in the 15 months following the eruption but then agree well with AVHRR to the end of 1993. An aspect of the initial eruption that is left out of the models is the injection of large volcanic ash particles, and this may partly explain the low AOD bias shortly after the eruption. While these large particles are short-lived (∼ 3 d; Guo et al., 2004a), their presence heats the volcanic cloud and can loft the volcanic plume to higher altitudes at rates of 1 km d⁻¹ in the first week, increasing the lifetime of the stratospheric sulfate and AOD (Stenchikov et al., 2021). The choice of injection may also be contributing to AOD underestimation, with E3SMv2-SPA and CESM2-WACCM both exhibiting lower tropical confinement than AVHRR and E3SMv2-presc (Fig. 3), leading to R_eff and lower AOD (Clyne et al., 2021). This emulates the behavior of other models with similar injection parameters (Quaglia et al., 2023). Slightly better tropical entrainment in E3SMv2-SPA is attributed to a higher plume compared to CESM2-WACCM (Fig. S6), so increasing the injection height or adding a lofting mechanism associated with volcanic ash may improve this feature in both E3SMv2-SPA and CESM2-WACCM.

The climate impacts of volcanic aerosol are explored through TOA flux (Fig. 4) and atmospheric temperature (Fig. 5) perturbations. In E3SMv2-SPA, the maximum reductions in all-sky ASR (−5.5 W m⁻²), OLR (−3.79 W m⁻²), and NET (−3.0 W m⁻²) flux compare reasonably well with those of ERBS (−5.3, −3.05, and −3.8 W m⁻², respectively). The clearest improvement of E3SMv2-SPA compared to ERBS is in ASR, which also has the most distinct Pinatubo signal in the data, especially during 1992 when aerosols are at their largest. There is also a stark difference between the E3SMv2-PA and E3SMv2-SPA, CESM2-WACCM, and E3SMv2-presc clear-sky and aerosol-only diagnostics, with an additional reduction of 2.5–3 W m⁻² in the three volcanic parameterized models. OLR is noisier but does show increased longwave absorption during the year 1992. Seasonal cloud albedo and surface cooling effects are introduced in models and observations through our choice of normalization year, and their effect can be removed in the models through comparison of clear-sky and aerosol-only TOA flux diagnostics. Temperature anomalies are calculated for the year 1992, corresponding to the greatest model spread in the TOA flux. Compared to MERRA-2 and RICH-obs, E3SMv2-SPA has a more comparable distribution of temperature anomalies than E3SMv2-PA, with notable improvements in the stratosphere over the midlatitudes. CESM2-WACCM is nearly identical to MERRA-2 temperature due to it being nudged to that particular dataset. This nudging puts CESM2-WACCM in very good agreement with MERRA-2, while E3SMv2-SPA overestimates 50 hPa temperatures and underestimates 850 hPa temperatures compared to MERRA-2. Aside from noting the improvement from E3SMv2-PA to E3SMv2-SPA, we caution against directly relating these temperature changes to aerosol size without a more rigorous removal of internal variability in the model (e.g., Santer et al., 2014).

The modeled aerosol size – represented by R_eff – is consistently underestimated compared to both in situ and remote sensing datasets (Fig. 6). E3SMv2-SPA also has slightly smaller R_eff values than CESM2-WACCM, where the slightly larger aerosols are attributed to the preferential growth of fewer initial Aitken-mode aerosols associated with the limited OH oxidation in CESM2-WACCM. Unpacking R_eff into its component modal distributions on a global scale illustrates the larger modal D_g and also smaller modal number concentration in CESM2-WACCM (Fig. 7). While useful for visualizing the evolution of a variety of volcanic eruptions over this period, validation with observations requires the finer-scale comparison of model and in situ data. The comparisons with balloon-borne WOPC data (Fig. 8) identify a key reason for the underestimation in R_eff, namely that modeled coarse-mode aerosols are not big enough in the models during the first year after Pinatubo. The aerosol modes reach a similar pattern of change after a year, and there is very close agreement in the microphysics of E3SMv2-SPA and CESM2-WACCM.

While the choice of bulk, sectional, or modal aerosol representation will impact R_eff (English et al., 2013; Laakso et al., 2022; Tilmes et al., 2023), this R_eff underestimation is present in modal as well as sectional models with similar injection parameters (e.g., ECHAM6-SALSA, SOCOL-AERv, UM-UKCA; Quaglia et al., 2023; Dhomse et al., 2020), indicating that it is not attributed solely to the choice of aerosol microphysical approach. We present a variety of potential reasons and solutions for why the model simulates aerosols that are too small, though testing these is beyond the scope of our study. As with AOD, higher injection levels and plume lofting could increase R_eff by increasing tropical confinement (Clyne et al., 2021). Improvement of the stratospheric dry bias in E3SMv2 may improve this disagreement by increasing aerosol size due to water uptake. There are also inherent model assumptions that may be contributing to the aerosol small bias: (1) the lack of interparticle van der Waals forces in both E3SMv2 and CESM-WACCM (though recent iterations of WACCM coupled with the CARMA sectional aerosol model do include these forces; Tilmes et al., 2023), which have been shown to drive aerosol nucleation and lead increased peak R_eff following Pinatubo-sized eruptions (English et al., 2013; Sekiya et al., 2016), and (2) the nucleation scheme used in these models possibly overestimating nucleation rates by 3–4 orders of magnitude (Yu et al., 2023), leading to more numerous smaller particles. The inability to create large enough aerosols in the models can affect the simulated climate impacts, and the size-based effects on scattering and absorption are explored through the single-particle Mie scattering.

The plotted Q_s (Figs. 9, 10) and Q_a (Fig. 13) tell an idealized story about how size affects the scattering of incoming solar radiation and absorption of outgoing terrestrial radiation. At the wavelength of peak incoming solar irradiance (0.5 µm; λ_solar), the larger observed aerosols manifest as weaker scatterers of incident radiation than E3SMv2-SPA and CESM2-WACCM. This suggests that the models will overestimate aerosol scattering, though the differences between model and observation are smaller for monthly as opposed to daily averages. For Q_a, the relationship is simpler, with a linear increase in aerosol absorption with increasing aerosol size. Our results show that the models may underestimate absorption of light at the wavelength of peak outgoing terrestrial radiation (10 µm; λ_earth) due to their smaller simulated aerosol size. When the effects of modeled interstitial aerosol water are included in the CESM2-WACCM refractive index, modeled Q_a is decreased further when compared to observations due to a lower absorption of λ_earth by water than by sulfate.

Care should be taken in interpreting Q_s and Q_a as there are some caveats to the above comparisons. While R_eff is a reasonable representation of aerosol size and is an important factor in aerosol optics, it is not the singular factor in aerosol optical properties. The aerosol number is neglected in these comparisons, which is a secondary determiner of AOD. The importance of the aerosol number can be seen in the diffuse radiation (Fig. 11), which has relative model behaviors that mirror those of AOD (Fig. 12). The number is also important in longwave absorption, where relative modeled longwave heating rates (Fig. 14) similarly mirror those in AOD. Another assumption in the calculation of scattering and absorption efficiencies is the assumption of a stratospheric sulfate refractive index. This is consistently done for remote sensing datasets where no additional information can be provided, but in models that prognostically calculate the volume-weighted soluble refractive index this may be misleading. Here, the same refractive index is assumed across observations and models to make them more comparable (Hess et al., 1998; Appendix B). It also avoids differences in CESM2-WACCM and E3SM that arise from the extremely low stratospheric water bias in E3SM. The accurate simulation of the number and refractive index are two aspects of this work that require further study.

It is unclear the extent to which biases in aerosol size and the resultant impacts on Q_a and Q_s are impacting the downstream climate impacts of stratospheric sulfate. While quantifying this impact and improving this bias are important next steps, we note that overall E3SMv2-SPA does well compared to observational datasets and produces Pinatubo sulfate burdens and AOD similar to CESM2-WACCM and other models run with comparable volcanic injection parameters (Quaglia et al., 2023). Furthermore, a free-running historical simulation run with E3SMv2-SPA recreates major historical volcanic AOD peaks compared to five-member E3SMv2 historical simulations with prescribed volcanic forcing (Fig. S1). While some disagreements exist, the general mirroring of the signals gives confidence in simulating variable injection parameters of varying magnitudes with this parameterization. Lastly, comparison of historical surface temperatures (Fig. S2) and net top-of-model radiative flux (Fig. S3) from two ensembles of a fully coupled E3SMv2-SPA historical simulation indicates a similar performance to E3SMv2 historical simulations with prescribed volcanic forcing. This indicates a reasonable climate representation in a historical setting.

This study has served to validate the stratospheric prognostic aerosol treatment in E3SMv2-SPA against observations while also showing it can have reasonable performance compared to the more comprehensive chemical treatment of CESM2-WACCM. While an interactive chemical treatment does have better overall agreement with observations of mass burden, AOD, and R_eff, E3SMv2-SPA can still serve as a viable alternative. We have connected the microphysical evolution of aerosols to radiative impacts, which then cascade into further downstream impacts (i.e., changes in temperature, precipitation, water vapor, etc.) that will be explored in follow-up papers. The use of a prescribed forcing dataset (E3SMv2-presc) does not allow the same connections to be drawn as there is no aerosol mass and size distribution that can evolve dynamically and feed back into the radiation balance in the Earth system, both directly and indirectly (i.e., aerosol–cloud interactions). Issues in representation and the lack of a dynamically evolving aerosol plume in prescribed forcing datasets could have a significant impact on detection and attribution studies looking at midlatitude agricultural productivity impacts, for instance, and illustrate the necessity of prognostically evolved aerosol for forthcoming attribution studies sensitive to spatio-temporal evolution. E3SMv2-SPA will facilitate prognostic aerosol simulations that include naturally produced stratospheric sulfate (i.e., volcanoes) as well as anthropogenic climate interventions to modify incoming solar radiation with stratospheric aerosol injection (SAI). Detection and attribution studies of societally relevant climate impacts in free-running fully coupled climate simulations with varied stratospheric sulfate aerosols are essential for casual identification of the source of impact and risk assessments (OSTP, 2023).

Appendix A: Effective radius and aerosol size distribution calculations

The aerosol effective radius (R_eff) – or area-weighted aerosol radius – is a good representation of a size distribution's optically relevant size (Hansen and Travis, 1974; Russell et al., 1996). This relies on the fact that intercepted light is proportional to the particle radius in Mie and geometric scattering regimes. These regimes can be designated by an aerosol size parameter, $x=\frac{\mathrm{2}\mathit{\pi }r}{\mathrm{wavlength}\left(\mathit{\lambda }\right)}$, with Mie scattering corresponding to $x\approx \mathrm{0.2}-\mathrm{2000}$ and geometric scattering corresponding to x>2000 (Petty, 2006). R_eff is calculated by taking the ratio of the third aerosol moment (volume) and second aerosol moment (cross-sectional area) (M₃ and M₂, respectively) of the aerosol size distribution (Hansen and Travis, 1974):

$$\begin{array}{}\text{(A1)}& R_{\mathrm{eff}}={\displaystyle \frac{M_{\mathrm{3}}}{M_{\mathrm{2}}}}.\end{array}$$

A lognormal size distribution is defined as

$$\begin{array}{}\text{(A2)}& n_N\left(\ln r\right)={\displaystyle \frac{\mathrm{d}N\left(r\right)}{\mathrm{d}\ln r}}={\displaystyle \frac{N}{\ln {\mathit{\sigma }}_{\mathrm{g}}\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{\left[-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\ln r-\ln r_{\mathrm{g}}}{\ln {\mathit{\sigma }}_{\mathrm{g}}}\right)}^{\mathrm{2}}\right]},\end{array}$$

where N is the total aerosol number concentration, r_g is the wet geometric mean radius, r is the radius, and σ_g is the geometric standard deviation. The analytic expression for the aerosol moments, M_n, of a lognormal distribution is

$$\begin{array}{}\text{(A3)}& M_n=N{r}_{\mathrm{g}}^n{e}^{\left[\frac{n^{\mathrm{2}}}{\mathrm{2}}{\left(\ln {\mathit{\sigma }}_{\mathrm{g}}\right)}^{\mathrm{2}}\right]}.\end{array}$$

For a multi-modal system, the modal moments (M_n,i) can be summed, giving the final equation for R_eff:

$$\begin{array}{}\text{(A4)}& R_{\mathrm{eff}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m{M}_{\mathrm{3},i}}{{\sum }_{i=\mathrm{1}}^m{M}_{\mathrm{2},i}}}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m{N}_i{r}_{g,i}^{\mathrm{3}}{e}^{\left[\frac{\mathrm{9}}{\mathrm{2}}\ln {{\mathit{\sigma }}_{g,i}}^{\mathrm{2}}\right]}}{{\sum }_{i=\mathrm{1}}^m{N}_i{r}_{g,i}^{\mathrm{2}}{e}^{\left[\mathrm{2}\ln {{\mathit{\sigma }}_{g,i}}^{\mathrm{2}}\right]}}}.\end{array}$$

When calculating R_eff over a range, z, of model hybrid or pressure levels (e.g., in the stratosphere, R_eff,strat), we integrate over the vertical levels weighted by the vertical layer thickness (h) and the aerosol surface area density (SAD).

$$\begin{array}{}\text{(A5)}& R_{\mathrm{eff},\mathrm{strat}}={\displaystyle \frac{{\sum }_z(R_{\mathrm{eff}}\times h\times \mathrm{SAD})}{{\sum }_z(h\times \mathrm{SAD})}}\end{array}$$

The thickness, h, is calculated using the hypsometric equation

$$\begin{array}{}\text{(A6)}& h={\displaystyle \frac{R_{\mathrm{d}}\times T}{g}}\ln {\displaystyle \frac{P_{z+\mathrm{1}}}{P_z}},\end{array}$$

and SAD (cm² cm⁻³) is represented by the second moment as

$$\begin{array}{}\text{(A7)}& \mathrm{SAD}=\mathrm{4}\mathit{\pi }{M}_{\mathrm{2}}.\end{array}$$

For stratospheric averages of lognormal aerosol distributions, we calculate the individual modal aerosol distributions (Eq. A2) and sum them across size bins (r) to create a multi-modal distribution at each grid cell (n_N(ln r)_3mode). We then weight this distribution by the grid cell thicknesses (h) and integrate above the tropopause:

$$\begin{array}{}\text{(A8)}& n_N{\left(\ln r\right)}_{\mathrm{strat}}={\displaystyle \frac{{\sum }_z(n_N{\left(\ln r\right)}_{\mathrm{3}\mathrm{mode}}\times h)}{{\sum }_z\left(h\right)}}.\end{array}$$

Appendix B: Single-particle Mie scattering

In order to understand how improvements to stratospheric aerosol size distributions affect simulated climate impact, aerosol scattering and absorption efficiencies (Q_s, Q_a) are calculated based on Mie theory and the solution of Maxwell's equations for the interaction of radiation with a sphere (Hansen and Travis, 1974). The Q_s (Q_a) is the surface area of the shadow cast by the particle due to scattering (absorption) of intercepted light divided by the geometric cross section (πr²). Q_s and Q_a are calculated using the Python package, miepython (https://miepython.readthedocs.io/en/latest/index.html#, last access: 20 June 2024). This package simulates absorption and scattering of incident light by non-absorbing or partially absorbing spheres and requires a complex aerosol refractive index ($n=a+bi$) and aerosol size parameter (aerosol circumference over the incident wavelength of light) as inputs. In the complex refractive index, the real (a) and imaginary (b) parts define scattering and absorption, respectively. Here, an effective size parameter is used as a proxy for the particle size distributions in our calculations:

$$\begin{array}{}\text{(B1)}& x_{\mathrm{eff}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }{R}_{\mathrm{eff}}}{\mathit{\lambda }}}.\end{array}$$

We note that because Mie theory calculates these values based on interactions with a single sphere and we are using the number-independent R_eff as an input, Q_s and Q_a do not account for effects of different number concentrations. The volume scattering and absorption coefficients (β_s, β_a; units of m⁻¹) do include the effect of aerosol number but are not generated based on our comparisons due to the ambiguity of assuming an effective aerosol number. Their equations are included below for completeness.

$$\begin{array}{}\text{(B2)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{s}}=N{Q}_{\mathrm{s}}\mathit{\pi }{r}^{\mathrm{2}}\text{(B3)}& {\displaystyle }{\mathit{\beta }}_{\mathrm{a}}=N{Q}_{\mathrm{a}}\mathit{\pi }{r}^{\mathrm{2}}\end{array}$$

The aerosol refractive index in the models is derived from a combination of sulfate aerosol and water due to the prevalence of sulfate in the stratosphere following Pinatubo. Because sulfate is a soluble species, a soluble refractive index (n_s) can be calculated with the volume mixing rule (Ghan and Zaveri, 2007):

$$\begin{array}{}\text{(B4)}& n_{\mathrm{s}}={\displaystyle \frac{\mathrm{1}}{V_{\mathrm{s}}}}\left({\displaystyle \frac{n_{\mathrm{wat}}{m}_{\mathrm{wat}}}{{\mathit{\rho }}_{\mathrm{wat}}}}+{\displaystyle \frac{n_{\mathrm{sulf}}{m}_{\mathrm{sulf}}}{{\mathit{\rho }}_{\mathrm{sulf}}}}\right),\end{array}$$

where n_wat, m_wat, and r_wat are the water complex refractive index, aerosol water mass (kg), and water density (kg m⁻³) and n_sulf, m_sulf, and r_sulf are the sulfate complex refractive index, sulfate mass (kg), and sulfate density (kg m⁻³). The m_sulf and m_wat are summed across the three aerosol modes. The soluble aerosol volume from sulfate and water across modes (V_s) is calculated as

$$\begin{array}{}\text{(B5)}& V_{\mathrm{s}}={\displaystyle \frac{m_{\mathrm{sulf}}}{{\mathit{\rho }}_{\mathrm{sulf}}}}+{\displaystyle \frac{m_{\mathrm{wat}}}{{\mathit{\rho }}_{\mathrm{wat}}}}.\end{array}$$

Sulfate and water refractive indices are taken from the literature (Hess et al., 1998; Hale and Querry, 1973). The sulfate refractive index reported in Hess et al. (1998) (n_Hess) is at 0 % relative humidity and corresponds to aerosol that is approximately 25 % water $/$ 75 % H₂SO₄ by mass (Seinfeld and Pandis, 2006). A similar solution ratio has been used in other remote sensing applications (see Sect. 3.1 and 3.5), and for our in situ and remote sensing datasets, we assume stratospheric aerosol optical properties are reasonably represented by n_Hess. The pure sulfate refractive index used in E3SM and CESM (n_sulf) is calculated from water and n_Hess following methods in Ghan and Zaveri (2007) and is nearly identical to n_Hess, albeit reported at a lower wavelength resolution. Given their similarity and the higher wavelength resolution in n_Hess, we use n_Hess instead of E3SM's pure sulfate refractive index for model single-particle scattering calculations.