Analysis of the global atmospheric background sulfur budget in a multi-model framework

Abstract

A growing number of general circulation models are adapting interactive sulfur and aerosol schemes to improve the representation of relevant physical and chemical processes and associated feedbacks. They are motivated by investigations of climate response to major volcanic eruptions and potential solar geoengineering scenarios. However, uncertainties in these schemes are not well constrained. Stratospheric sulfate is modulated by emissions of sulfur-containing species of anthropogenic and natural origin, including volcanic activity. While the effects of volcanic eruptions have been studied in the framework of global model intercomparisons, the background conditions of the sulfur cycle have not been addressed in such a way. Here, we fill this gap by analyzing the distribution of the main sulfur species in nine global atmospheric aerosol models for a volcanically quiescent period. We use observational data to evaluate model results. Overall, models agree that the three dominant sulfur species in terms of burdens (sulfate aerosol, OCS, and SO${}_{\mathrm{2}}$) make up about 98 $\mathit{\%}$ stratospheric sulfur and 95 $\mathit{\%}$ tropospheric sulfur. However, models vary considerably in the partitioning between these species. Models agree that anthropogenic emission of SO${}_{\mathrm{2}}$ strongly affects the sulfate aerosol burden in the northern hemispheric troposphere, while its importance is very uncertain in other regions, where emissions are much lower. Sulfate aerosol is the main deposited species in all models, but the values deviate by a factor of 2. Additionally, the partitioning between wet and dry deposition fluxes is highly model dependent. Inter-model variability in the sulfur species is low in the tropics and increases towards the poles. Differences are largest in the dynamically active northern hemispheric extratropical region and could be attributed to the representation of the stratospheric circulation. The differences in the atmospheric sulfur budget among the models arise from the representation of both chemical and dynamical processes, whose interplay complicates the bias attribution. Several problematic points identified for individual models are related to the specifics of the chemistry schemes, model resolution, and representation of cross-tropopause transport in the extratropics. Further model intercomparison research is needed with a focus on the clarification of the reasons for biases, given the importance of this topic for the stratospheric aerosol injection studies.

1 Introduction

Sulfur in the atmosphere modulates incoming solar radiation, affects the ozone layer, fertilizes soils, and impacts air quality in industrial areas. The most abundant gaseous sulfur species in the atmosphere are carbonyl sulfide (OCS) and sulfur dioxide (SO${}_{\mathrm{2}}$). Shorter-lived or emitted in smaller amounts, and therefore less abundant, are dimethyl sulfide (DMS) emitted from marine phytoplankton, hydrogen sulfide (H${}_{\mathrm{2}}$S), or carbon disulfide (CS${}_{\mathrm{2}}$) . Only a fraction of the sulfur emitted at the surface is transported to the stratosphere, with the majority scavenged in the mid-troposphere e.g.,. In the stratosphere, these sulfur-containing species get photolyzed and oxidized to eventually form sulfuric acid (H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$), the final oxidation product. Because of its low saturation vapor pressure, gaseous H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ then readily condenses and/or nucleates in combination with water vapor to aerosol particles, forming the “Junge layer”, a layer of aqueous sulfuric acid droplets (in short, “sulfate aerosol”) in the region between the tropopause and about 10 hPa . During volcanically quiescent (background) periods, the Junge layer is maintained by surface emissions of these precursor gases and their oxidation products and is assumed to be relatively constant. On the other hand, with the injection of wildfire smoke and the influence of frequent small and moderate volcanic eruptions, there are only a few years within the satellite era, when the stratospheric aerosol layer can be considered close to background or unperturbed e.g.,.

Most of the research related to the aerosol layer has been focused on large volcanic eruptions and their influence on climate e.g.,, atmospheric composition e.g.,, and dynamics e.g.,. While large volcanic events are one of the main natural climate drivers , small volcanic eruptions have also been shown to significantly contribute to the global radiative forcing and climate variability . In addition, the background aerosol layer itself undergoes substantial inter-annual variations . The stratospheric aerosol layer has become of interest for a more controversial reason as well; to moderate global climate warming, it has been proposed to inject sulfate aerosol precursors in the stratosphere in an attempt to mimic the global surface cooling generated by large volcanic eruptions and thus counteract the climate warming from increased greenhouse gases e.g.,. Predicting the effects of the stratospheric aerosol variations requires the simulation of multiple coupled processes with complex global general circulation models (GCMs) that are still subject to significant uncertainties. In these global models, chemical species and aerosols can be either prescribed or calculated interactively as prognostic variables. The former approach, while being computationally less expensive, is limited by uncertainties in the observations used to derive the prescribed distributions and does not account for the coupling of processes and internal feedbacks that would impact the distributions themselves. Furthermore, biases of up to 20 % in aerosol extinction measurements across different satellite instruments mean that small variations cannot currently be adequately quantified by observations . Models with interactive aerosol schemes and chemistry, on the other hand, have many parameters, and potentially more degrees of freedom (and therefore more sources of uncertainty) but can account explicitly for the feedbacks between aerosol microphysics and dynamical and chemical processes.

With the growing availability of computational resources and scientific evidence of a potentially large role of the Junge layer in future climate , increasing numbers of GCMs are now including interactive sulfur and aerosol schemes to improve the representation of relevant chemical processes and associated feedbacks. To evaluate the individual model performances and characterize the inter-model uncertainty in the involved processes, there have been several model intercomparison studies focused on elevated aerosol conditions due to volcanic events and artificial sulfur injections . So far, apart from the limited and quite old global aerosol model intercomparison for non-volcanic conditions described in the report, all previous global stratospheric aerosol model intercomparison studies have focused on the volcanically perturbed aerosol layer. This is an unusual situation because, whatever processes or climate components are considered in a model evaluation, models are generally assessed first for background conditions before moving to perturbed conditions. The few model studies on the background state of the aerosol layer are almost all single model studies, leaving the possibility that some of the results and conclusions might be model-dependent. The results of these studies show quite good agreement with observations for specific parameters but also reveal discrepancies for others . In the majority of models, the background conditions have not been evaluated at all, and a comprehensive and extensive multi-model assessment with interactive chemistry schemes for all sulfate aerosol precursor gases in the background state is still pending. Such a study has the potential to reveal common deficiencies in the model representation of specific processes, especially concerning differences between the models, which are hard to identify under volcanically perturbed conditions but still have repercussions for model performance. For example, noted a difference in the aerosol effective radius among models in experiments on the 1991 Pinatubo eruption, which cannot be addressed in detail in a perturbed state. Furthermore, showed, with measurements from the Stratospheric Aerosol and Gas Experiment on the International Space Station (SAGE III/ISS; ), that small volcanic perturbations of the background aerosol layer can lead to an increase or reduction in the aerosol effective radius, depending on the regional background conditions of the individual events. Finally, characterizing the background state and its modeling uncertainties can be useful for the next Coupled Model Intercomparison Project (CMIP) phase preparations, as the semi-background aerosol state (averaged 1850–2014) is usually used for the Diagnostic, Evaluation and Characterization of Klima (DECK) and Shared Socioeconomic Pathway (SSP) experiments .

The background (BG) aerosol in the stratosphere is highly dependent not only on the precursor gases and background chemistry e.g., but also on the variability and evolution of atmospheric dynamics, which controls the stratosphere–troposphere exchange, as well as the general stratospheric circulation, the so-called Brewer–Dobson circulation (BDC) . Thus, the model performance in terms of background aerosol layer climatology and variability can be expected to be affected by underlying model transport biases. This sensitivity of the modeled aerosol layer to the stratospheric transport is more difficult to assess for volcanically perturbed conditions because sulfate aerosols are much larger than for background conditions, and hence, sedimentation plays a larger role in terms of aerosol transport and global redistribution. The issue of dynamical differences between models has been highlighted in several studies. For example, show large inter-model differences in mixing activity (including horizontal and vertical mixing, as well as vertical diffusion), which affects the age of air (AoA) and therefore stratospheric transport of chemical species. A recent Coupled Model Intercomparison Project phase 6 (CMIP6) model evaluation of the BDC has revealed that, while models generally agree on the AoA in the lower branch, larger differences exist in the middle and upper stratosphere . also show how a coarser model resolution negatively impacts the representation of the tropical and polar vortex transport barriers. Similarly, show that increasing the model vertical resolution strengthens the sub-tropical transport barrier, increasing the residence time of chemical species or aerosol in the tropics. show that the stratospheric aerosol layer is also highly modulated by the quasi-biennial oscillation (QBO), with non-linear QBO-phase-dependent regional effects.

To date, no extensive model intercomparison exists on background atmospheric sulfur burdens and distributions and characterization of related uncertainties. This activity has been proposed by in the framework of the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA–MIP), among other experiments, to comprehensively study and intercompare the representation of stratospheric aerosol processes in different models. The ISA–MIP BG experiment is designed to reproduce volcanically unperturbed conditions of the Junge layer in a 20-year time slice simulation with predefined boundary conditions. This experiment is expected to reveal common model deficiencies that are not visible in volcanic experiments, providing valuable information for guiding improvements in stratospheric aerosol models. Here, we follow the proposed BG setup and compare simulations from nine atmospheric models participating in ISA–MIP, as listed in Table . The aim is to quantify the range of simulated burdens and distributions of stratospheric aerosols and to evaluate the model results against satellite-derived observations. Furthermore, we identify existing uncertainties in state-of-the-art stratospheric aerosol models.

The experimental setup and the models involved are described in Sect.  and , respectively. Section  describes all observational datasets used for model evaluation. In Sect. , we discuss the full sulfur budget in the participating models, including some of the main chemical processes influencing the aerosol layer, emissions of aerosol precursors, cross-tropopause fluxes, stratospheric burdens, reactions, and deposition. We also present the total sulfur budget summed over all sulfur species. Seasonal cycles and meridional distributions are presented in Sect. . The distribution of three major sulfur species is shown in Sect. . Section  discusses effective radius and surface area density. Finally, the conclusions of this work are presented in Sect. .

2 Data and methods

2.1 Experimental setup

We follow the setup described by in the BG experiment BG_QBO for all models, unless otherwise specified in the model description (hereafter called REF). The simulations are set up as 20-year time slice simulations using repeating boundary conditions of the year 2000. A second simulation, termed BG_NAT (hereafter called NAT), has the same setup, except that all anthropogenic sulfur emissions were excluded . The difference in SO${}_{\mathrm{2}}$ emissions between these two experiments is shown in Fig. . All aerosol and sources of aerosol precursors, except explosive volcanic eruptions, are included in this study. In the models without an internally generated QBO (see Table ), the QBO is nudged to the 1981–2000 period. All models have been recommended to use prescribed sea surface temperatures (SSTs) and sea ice coverage (SIC) from the Met Office Hadley Center Observational Dataset . Sulfur emissions from anthropogenic sources and biomass burning are taken from the Monitoring Atmospheric Composition and Climate (MACCity) inventory , repeating the year 2000. Emissions from continuously degassing volcanoes are given by , based on . OCS surface concentrations are prescribed and constant at 510 pptv (parts per trillion by volume) . If allowed by the model setup, DMS emissions are calculated online using concentrations in the global oceans given by . Some models include similar online calculations for dust and sea salt, in which case, the oceanic concentrations for these compounds are also taken from . Otherwise, models use their usual database. Stratospheric burdens, if not directly provided in the output (in CAM5–CARMA, ECHAM5–HAM, WACCM6–CARMA, WACCM6–MAM4, and UM–UKCA), were calculated from monthly mean mixing ratios, using standard air in CAM5–CARMA and UM–UKCA, and the provided air mass in all other models. The stratospheric burden was then calculated by masking out all grid boxes below the model tropopause, not accounting for the volume of partially stratospheric grid boxes.

2.2 ISA–MIP models

2.2.1 SOCOL–AERv2

The global atmosphere-only chemistry–climate model SOCOL–AERv2 consists of the interactively coupled dynamical core ECHAM5 and the chemistry model MEZON, forming SOCOLv3 , as well as the aerosol model AER . The model uses a triangular truncation at wavenumber 42 (T42), which corresponds to a resolution of about 2.8° $\times$ 2.8° and it extends vertically to 0.01 hPa (or about 80 km). In this study, we use a vertical resolution of 39 levels. SOCOL–AERv2 uses a sectional aerosol scheme, differentiating 40 size bins. All aerosol in this model is pure sulfate aerosol. The microphysics in SOCOL–AERv2 include a nucleation scheme by , condensation and evaporation according to and , and coagulation . Sedimentation occurs according to and , whereas aerosol composition is derived from . SOCOL–AERv2 uses interactive dry and wet deposition schemes based on the DRYDEP and SCAV modules. In total, SOCOL–AERv2 distinguishes eight sulfur species (OCS, CS${}_{\mathrm{2}}$, MSA, DMS, H${}_{\mathrm{2}}$S, SO${}_{\mathrm{2}}$, SO${}_{\mathrm{3}}$, and H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$), as well as sulfate aerosol and 27 reactions, including sulfur species. In addition to sulfur species, the model also includes oxygen, hydrogen, nitrogen, carbon, chlorine, and bromine species .

Table 1

A list of models that participated in this study and are part of the ISA–MIP${}^{*}$ project. Described are the horizontal and vertical resolutions, as well as if the QBO is internally generated. This is followed by a short description of the aerosol scheme. CAM5–CARMA and WACCM5–CARMA have two sets of aerosol bins for different aerosol types (denoted by 2 $\times$ 20).

Model	Lat	Long	Levels	Top	QBO	Aerosol scheme	References
						(modes/bins)
SOCOL–AERv2	2.8$\mathrm{°}$	2.8$\mathrm{°}$	39	80 km	Nudged	sectional (40)	,
ULAQ–CCM	5$\mathrm{°}$	6$\mathrm{°}$	126	70 km	Nudged	sectional (22)	,
ECHAM6–SALSA	1.9$\mathrm{°}$	1.9$\mathrm{°}$	95	80 km	Internally gen.	sectional ($\mathrm{10}+\mathrm{7}$)
CAM5–CARMA	1.9$\mathrm{°}$	2.5$\mathrm{°}$	56	40 km	Nudged	sectional ($\mathrm{2}\times \mathrm{20}$)
ECHAM5–HAM	2.8$\mathrm{°}$	2.8$\mathrm{°}$	90	80 km	Internally gen.	modal (7)
WACCM6–CARMA	1.9$\mathrm{°}$	2.5$\mathrm{°}$	70	140 km	Nudged	sectional ($\mathrm{2}\times \mathrm{20}$)
WACCM6–MAM4	1.9$\mathrm{°}$	2.5$\mathrm{°}$	70	140 km	Nudged	modal (4)
MIROC–CHASER	2.8$\mathrm{°}$	2.8$\mathrm{°}$	57	52 km	Nudged	modal (3)	,
UM–UKCA	1.875$\mathrm{°}$	1.25$\mathrm{°}$	85	85 km	Internally gen.	modal (7)	,
							,

${}^{*}$ https://isamip.eu/home (last access: 23 April 2024)

2.2.2

ULAQ–CCM

The global-scale chemistry–climate model ULAQ–CCM (University of L'Aquila Chemistry Climate Model) has a resolution of 5° $\times$ 6° (T21) and uses 126 log-pressure levels, reaching from the Earth's surface to 0.04 hPa. It treats sulfate, organic and black carbon, dust, sea salt, nitrate, and polar stratospheric cloud (PSC) aerosols . Each type of aerosol is treated separately in terms of surface fluxes, transport, and removal from the atmosphere. The wet and dry deposition schemes are based on . Included in the chemistry module are species from the O${}_x$, NO${}_y$, NO${}_x$, CHO${}_x$, HO${}_x$, Cl${}_y$, Br${}_y$, and SO${}_x$ families. This includes six sulfur species: OCS, CS${}_{\mathrm{2}}$, DMS, H${}_{\mathrm{2}}$S, H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, and SO${}_{\mathrm{2}}$, as well as the long-lived species N${}_{\mathrm{2}}$O, CH${}_{\mathrm{4}}$, CO, hydrocarbons, CFCs, hydrochlorofluorocarbons, and halons. . For ULAQ, only 10 years of the simulation were conducted.

2.2.3

WACCM6–CARMA and WACCM6–MAM4

In this study, we include simulations with the Community Earth System Model version 2 in its high-top configuration, named the Whole Atmosphere Community Climate Model version 6 (WACCM6). We use the middle atmosphere (MA) chemistry mechanism . Hereafter, we call this model setup WACCM6–MA. In our setup, WACCM6–MA has 70 vertical levels reaching up to 140 km above the surface. We set the horizontal resolution of 1.9$\mathrm{°}$ $\times$ 2.5$\mathrm{°}$. This model includes a comprehensive chemistry scheme in the stratosphere, mesosphere, and lower thermosphere, while only representing limited chemistry in the troposphere . The model accounts for sulfur chemistry of important precursor emissions for both the troposphere and stratosphere, including four sulfur species, namely OCS, DMS, SO${}_{\mathrm{2}}$, and H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ .

WACCM6 is coupled to two different aerosol microphysical modules. The Modal Aerosol Model (MAM4) , which includes updated prognostic stratospheric sulfate aerosols , is the default aerosol scheme used in CAM-chem and WACCM6 of CESM2. Four modes are described by MAM4 microphysics: Aitken, accumulation, and coarse modes, as well as a primary carbon mode . The geometric standard deviation of MAM4 for the Aitken and accumulation mode is 1.6, while for the coarse mode, it is 1.2 .

The second aerosol microphysical model coupled to WACCM6 is the Community Aerosol and Radiation Model for Atmospheres (CARMA) version 4.0, which enables size-resolved or sectional cloud droplets and aerosol particles . The CARMA aerosol model includes prognostic aerosols for both the troposphere and the stratosphere, as discussed in , and additional changes are highlighted in . Apart from pure sulfate aerosol, WACCM6–CARMA is one of two models participating in this study which includes an internally mixed group of aerosol. It involves sulfate, primary and secondary organics, black carbon, dust, and sea salt. The model divides each group into 20 discrete mass bins, as defined by . The mixed aerosol group specifies bins with radii between 0.05 and 8.7 $\mathrm{\mu }\mathrm{m}$, whereas radii of the pure sulfate group range from 0.2 nm to 1.3 $\mathrm{\mu }\mathrm{m}$. The aerosol composition is based on .

2.2.4

ECHAM6–SALSA

The aerosol–climate model ECHAM6–SALSA (ECHAM6.3–HAM2.3–MOZ1.0–SALSA2.0 is comprised of the ECHAM6.3 general circulation model and the HAM aerosol module . The last component is the aerosol microphysics module SALSA2.0 . The model was set up with a T63 resolution, corresponding to a $\mathrm{1.9}\times \mathrm{1.9}$ horizontal grid. Furthermore, it uses 95 vertical levels with a top at 0.01 hPa. The microphysical scheme SALSA uses 10 fixed size bins, ranging from 3 nm to 10 $\mathrm{\mu }\mathrm{m}$, while the seven largest bins additionally treat soluble and insoluble aerosol . In this study, we use the parameterized sulfuric acid–water binary homogeneous nucleation parameterization for nucleation. The analytical predictor of condensation (APC) scheme is applied to calculate condensation , while coagulation is treated according to . Apart from sulfate, SALSA also includes organic aerosol, sea salt, black carbon, and mineral dust. The deposition and sedimentation in SALSA are presented by . ECHAM6–SALSA uses a simplified chemistry scheme from HAM and includes the oxidation of DMS and SO${}_{\mathrm{2}}$ via a range of oxidizing agents (OH, H${}_{\mathrm{2}}$O${}_{\mathrm{2}}$, NO${}_{\mathrm{2}}$, and O${}_{\mathrm{3}}$) prescribed by a monthly mean climatology. ECHAM6–SALSA includes three of the main sulfur gases (DMS, SO${}_{\mathrm{2}}$, and H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$), whereas OCS is not included.

2.2.5

ECHAM5–HAM

ECHAM5–HAM uses the high-top version of ECHAM5 and is coupled to HAM, an aerosol microphysical model . The horizontal grid has a 2.8° $\times$ 2.8° resolution, whereas vertically, there are 90 layers up to 0.01 hPa, corresponding to 80 km. Microphysics in HAM treats the oxidation of sulfur, including sulfate aerosol formation. This encompasses nucleation, accumulation, condensation/evaporation, and coagulation. To improve the stratospheric aerosol representation, modifications were made to the microphysical core M7 of HAM , especially for high sulfur loads after volcanic eruptions. HAM uses a modal size distribution comprised of four modes. The simulations for this paper used the nucleation, Aitken, and accumulation mode with a mode width of 1.59 and a coarse mode with a mode width of 2. Another addition was made to HAM in the form of a simple stratospheric sulfur chemistry scheme . The sulfur chemistry in ECHAM5–HAM tracks four sulfur gases, namely OCS, DMS, SO${}_{\mathrm{2}}$, and H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$. As the chemistry scheme is not fully interactive, monthly fields for OH, NO${}_{\mathrm{2}}$, and O${}_{\mathrm{3}}$, as well as photolysis rates of OCS, H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, SO${}_{\mathrm{2}}$, and O${}_{\mathrm{3}}$ are prescribed on a monthly mean basis. A general description of the performance of HAM is described in , , and .

2.2.6 CAM5–CARMA

CAM5–CARMA is a low-top version of the Community Earth System Model (CESM1), coupled to the aerosol microphysical model CARMA, which is described in Sect. . It has a horizontal resolution of 1.9° $\times$ 2.5° and runs on hybrid 56 vertical levels . The model includes full stratospheric and tropospheric chemistry, using the chemistry module MOZART-4 . CARMA tracks organic carbon, black carbon, dust, and sea salt, as well as an internally mixed type . Secondary organic aerosol is included and based on . CARMA provides a sectional aerosol scheme, tracking 20 particle size bins for aerosol and another 20 for mixed aerosol. DMS emissions are based on . The chemistry scheme in CAM5–CARMA includes 230 chemical reactions. Sulfur chemistry is based on and includes 22 gas-phase and 5 heterogeneous reactions, summarized in . Three sulfur species are tracked in CAM5–CARMA, additionally to sulfate aerosol, namely OCS, SO${}_{\mathrm{2}}$, and H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$. For this simulation, CAM5–CARMA was nudged to the MERRA-2 reanalysis. Instead of the 20-year time slice simulation, we use 20 ensemble members for the year 2000.

2.2.7 MIROC–CHASER

The global chemistry–climate model MIROC–CHASER consists of the Model for Interdisciplinary Research on Climate (MIROC) and the atmospheric chemistry model CHASER and the Spectral Radiation–Transport Model for Aerosol Species (SPRINTARS) . For this study, the model is set up with a 2.8° $\times$ 2.8° horizontal resolution and 57 vertical levels up to 52 km. The aerosol module SPRINTARS tracks sulfate aerosol with three modes (nucleation, Aitken, and accumulation) and uses the bulk approach for black carbon and organic matter, dust, and sea salt . DMS emissions in MIROC–CHASER are a function of downwelling short-wave radiation. Nucleation is based on , while coagulation follows the same scheme as ECHAM5–HAM . The chemical scheme in CHASER includes 93 species, as well as 263 reactions. Sulfur chemistry is included in the form of 12 reactions, as well as 4 of the main sulfur species (SO${}_{\mathrm{2}}$, SO${}_{\mathrm{4}}$, DMS, OCS) .

2.2.8 GA4-UM–UKCA

The first GA4-UM–UKCA (hereafter UM–UKCA) simulation submitted for the ISA–MIP BG$\mathit{\_}$QBO (here called REF) experiment uses the interactive stratospheric aerosol configuration of the UM–UKCA model . The model runs with a horizontal resolution of 1.875° $\times$ 1.25° and on 85 levels, with a model top at approximately 85 km. Specifically, this first UM–UKCA submission to BG applies the identical version 3 (v3) stratosphere–troposphere UKCA codebase also run for the ISA–MIP HErSEA–Pinatubo experiment ; as further analyzed by , each of the ISA–MIP UM–UKCA runs within the GA4 configuration of the UM general circulation model .

This version 3 of stratosphere–troposphere UM–UKCA comprises version 8.2 of the GLOMAP-mode aerosol microphysics module see implemented within the RJ4.0 configuration of the UK Chemistry and Aerosol sub-model, as released to the UK academic community within GA4 . The chemical scheme accounts for OCS, SO${}_{\mathrm{2}}$, DMS, H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, SO${}_{\mathrm{3}}$, and sulfate aerosol . For the REF experiment, the simulations are within the year 2000 time slice atmosphere-only simulations, with boundary conditions and tropospheric chemistry and aerosol emissions identical to that described by and with the aerosol–radiation and aerosol–cloud interaction radiative effects and UM–UKCA simulated tropospheric and stratospheric ozone layers fully interactive with the radiative transfer module within GA4.

The 20-year simulation analyzed is the last 5 years from the UM–UKCA v3 simulation shown in , with an extension for a further 15 years for the REF experiment. As explained in , v3 UM–UKCA includes heterogeneously nucleated sulfuric acid aerosol particles with the 7.9 t d${}^{-\mathrm{1}}$ meteoric smoke particle (MSP) climatology (v3 low MSP). As configured for the analysis in , the simulation has the modal desert dust emissions switched off, with desert dust radiative effects from the UM sectional interactive dust scheme rather than from GLOMAP-mode.

The TS2000 atmosphere-only RJ4.0 UM–UKCA model used here is identical to that also applied for the 2000 volcanic forcing perturbed-parameter ensemble and equivalent also to the pre-industrial setting GA4 UM–UKCA v3 simulations for the Volcanic Forcings Model Intercomparison Project (VolMIP) interactive Tambora experiment .

2.3 Observational and reanalysis data

We use several satellite-derived datasets for the stratospheric aerosol layer , stratospheric SO${}_{\mathrm{2}}$ , stratospheric OCS , and in situ measurements from high-altitude balloon soundings . An overview of these datasets is provided in Table .

For each of these datasets, we use the period 1999–2004 (unless otherwise specified) that is relatively unaffected by volcanic eruptions. Comparing transient observations with time slice simulations comes with certain caveats, including the influence of the QBO. The influence of the QBO is not limited to the tropical stratosphere but also affects export from the tropics to mid-latitudes and modulates the strength of the winter polar vortex . Therefore, we do not try to validate models against observational data. Rather, we attempt to qualitatively assess the differences and refer the reader to previous papers in which quantitative model validations have been conducted.

Table 2

Observational and reanalysis datasets used in this study and the original resolution of each dataset.

Variable	Dataset/instrument	Lat	Long	Levels	Extent	References
Aerosol burden	SAGE-3$\mathit{\lambda }$	36°	1°	70	5–40 km	,
						,
						Beiping Luo (personal communication, 2023)
Optical depth and extinction	GloSSAC V2.2	32°	1°	80	0–40 km
SO${}_{\mathrm{2}}$	MIPAS	18°	6°	36${}^{*}$	10–45 km
OCS	MIPAS	18°	6°	26${}^{*}$	10–35 km
Surface area density	SAGE II	32°	1°	–/–	0–40 km	,
and effective radius
Temperature and zonal winds	ERA-Interim	64°	128°	60	0–0.1 hPa
Effective radius	OPC Laramie	41° N	105° W	–/–	$\sim$ 13–30 km

${}^{*}$ The numbers of levels for MIPAS do not reflect the vertical resolution. These are provided in the dataset and the related publications.

We use an aerosol dataset derived using the 3-$\mathit{\lambda }$ mechanism, as described by , which is based on the Global Space-based Stratospheric Aerosol Climatology version 2.2 (GloSSACv2.2) dataset (see below). These data are derived from extinction measurements at three wavelengths by the Stratospheric Aerosol and Gas Experiment II (SAGE II; ). The dataset, hereafter called SAGE-3$\mathit{\lambda }$, includes the zonal mean distribution of aqueous H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ concentrations as monthly mean values from 1979–2021 at 36 latitudes, where the data beyond 80° N and S have been extrapolated.

Beiping Luo, personal communication, 2023.

SAGE-3$\mathit{\lambda }$ provides data at altitudes every 500 m between 5 km and 40 km. We apply this dataset above the lapse rate tropopause derived from ERA-Interim (ERA-I) temperature data .

We evaluate the model-derived extinctions at $\sim$ 525 nm (the wavelength corresponding to SAGE II data with the smallest uncertainty) and compare them with the GloSSACv2.2 data, although optical properties are not the main focus of this analysis. The GloSSACv2.2 dataset includes the solar occultation SAGE II measurements at four wavelengths out of the seven SAGE II wavelengths available . The dataset is composed of monthly zonal mean values at 32 latitudes, ranging from 80° N to 80° S and vertical levels from 0.5 to 40 km at 500 m intervals.

A dedicated level-3 dataset for SO${}_{\mathrm{2}}$, derived from the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS), is used in the present study. To achieve global distributions of SO${}_{\mathrm{2}}$ with a vertical coverage from 10 to 45 km we have combined two datasets: (1) the gridded dataset used in , based on the MIPAS single retrievals from 10 to 23 km altitude and 5 d binning, and (2) the MIPAS monthly mean retrievals, with a vertical coverage of 15–45 km . To minimize the volcanic influence, several months in the dataset that appeared to be affected by eruptions were excluded from the analysis, as described by . As this leaves significant gaps in the dataset, we use an average over the whole period from 2002 to 2012 for our analysis. Furthermore, we use the temperature distribution from ERA-I to derive the lapse rate tropopause and calculate stratospheric burdens.

For OCS, we utilized a gridded dataset based on the MIPAS retrievals by with a three-dimensional sampling of 10° latitude, 60° longitude, and 1 km altitude with a vertical coverage of 10 to 35 km and a temporal averaging of 5 d. The temporal coverage of this dataset from July 2002 to April 2012 does not include the first 3.5 years of our chosen time period. However, since OCS does not significantly depend on volcanic emissions and trends in background OCS are small, these influences are disregarded here . To calculate stratospheric burdens, the ERA-I-derived tropopause was applied.

We also use the effective radius and aerosol surface area density dataset derived from SAGE II with the SAGE version 7.0 algorithm . Aerosol was assumed to be composed of aqueous sulfuric acid solution at 75 wt $\mathit{\%}$ H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$.

We also compare the simulated size distributions to balloon measurements from Laramie, Wyoming, at 41° N and 104° W . The data have been collected using optical particle counters (OPCs), which measure the particle concentration in 12 size bins from 0.15–2 $\mathrm{\mu }\mathrm{m}$ .

3 Results and discussion

3.1 The global atmospheric sulfur budget

In Fig. , we present the global atmospheric sulfur budget represented by burdens and fluxes of major sulfur species, based on the input by the nine global atmospheric aerosol models involved in this study (see Table ). The major sulfur (S) species are specified in gigagrams of sulfur of each sulfur species ($\mathrm{Gg}\left(\mathrm{S}\right)$) for the burdens and gigagrams of sulfur of each sulfur species per year ($\mathrm{Gg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$) for fluxes. For the four models CAM5–CARMA, WACCM6–CARMA, WACCM6–MAM4, and UM–UKCA, stratospheric burdens are calculated from monthly mean volume or mass mixing ratios integrated from the model tropopause upward, as these burdens were not directly provided by some of the modeling groups. Mean emission and deposition fluxes are calculated by averaging the corresponding output from those models that provided them.

Figure 1

The atmospheric sulfur cycle with burdens (in $\mathrm{Gg}\left(\mathrm{S}\right)$) and fluxes ($\mathrm{Gg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$) for S gases and sulfate aerosols (displayed in analogy to , and ). The term “aerosol” refers only to sulfate aerosol. The figure includes all burdens as provided by the models involved. For CAM5–CARMA, WACCM6–CARMA, WACCM6–MAM4, and UM–UKCA, stratospheric burdens are calculated from the monthly mean volume or mass mixing ratios integrated from the model tropopause upward because these burden values were not provided. Burdens and fluxes are averaged over the entire period. We use SAGE-3$\mathit{\lambda }$ for aerosol and MIPAS for SO${}_{\mathrm{2}}$ and DMS emissions from , which are all observation-derived. Additionally, we show the aerosol wet depositions from a multi-model mean from the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP) , as previously used in . These observation-derived data and the ACCMIP values are in parentheses. Most fluxes are given as one-way fluxes, while cross-tropopause transport values, as well as OCS and H${}_{\mathrm{2}}$S emissions, are provided as net fluxes. Fluxes are only available for one or two models at a time (SOCOL–AERv2 and ULAQ–CCM or SOCOL–AERv2 and ECHAM5–HAM). Emissions of OCS and H${}_{\mathrm{2}}$S were prescribed as mixing ratios in the surface layer in some models. In that case, emission and deposition net fluxes for both species are derived by balancing the sum of the other fluxes of all species.

[Figure omitted. See PDF]

Table 3

The multi-model mean, standard deviation, and coefficient of variation in all major species in the troposphere, the stratosphere, and their emissions and depositions. The total tropospheric burden, total stratospheric burden, total emission, and total deposition include the values from all species and the fraction of each compound of the total mass of sulfur in each category.

	Tropospheric burden ($\mathrm{Gg}\left(\mathrm{S}\right)$)							Stratospheric burden ($\mathrm{Gg}\left(\mathrm{S}\right)$)							Emission ($\mathrm{Tg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$)					Deposition ($\mathrm{Tg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$)
	OCS	DMS	SO${}_{\mathrm{2}}$	H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$	sulfate	total		OCS	DMS	SO${}_{\mathrm{2}}$	H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$	sulfate	total		SO${}_{\mathrm{2}}$	DMS	other${}^{\mathrm{a}}$	total		SO${}_{\mathrm{2}}$	sulfate	other${}^{\mathrm{b}}$	total
Model mean	2391	99	347	0.9	613	3364		319	0.06	12	1.9	156	486		65	28	4.3	98		33	61	3.9	98
Standard deviation	100	69	123	0.83	161	282		50	0.05	6.8	1.8	51	85		4	4		2		14	21		18
Coeff. of variation	4 %	70 %	35 %	95 %	26 %	8 %		16 %	82 %	56 %	92 %	33 %	18 %		7 %	16 %		2 %		48 %	44 %		21 %
Fraction of total S	68 %	3 %	10 %	0.02 %	17 %			64 %	0.01 %	2.5 %	0.4 %	31 %			65 %	28 %				37 %	61 %

${}^{\mathrm{a}}$ Available only from SOCOL–AERv2 by summing up the chemical fluxes from H${}_{\mathrm{2}}$S, CS${}_{\mathrm{2}}$, and OCS to SO${}_{\mathrm{2}}$. ${}^{\mathrm{b}}$ DMS, H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, and MSA, with DMS being based on MIROC–CHASER and WACCM6–CARMA and MSA being based only on SOCOL–AERv2.

Observational data, as well as one ACCMIP estimate (in parentheses), are from the SAGE-3$\mathit{\lambda }$ and MIPAS datasets described in Sect.  for stratospheric sulfate aerosol, SO${}_{\mathrm{2}}$, OCS, and DMS emissions, respectively, and the aerosol wet deposition rates. Figure  is constructed in analogy to the figures of the sulfur cycle shown by and . Subsequently, we refer to sulfate aerosol often simply as aerosol. Additionally, we present the multi-model mean, standard deviation, and coefficient of variation in Table . The coefficients of variation are calculated as the standard deviation divided by the model mean. Mean burdens are calculated by averaging the output of all nine models (except OCS, which is an average over eight models without ECHAM6–SALSA, which does not treat OCS).

3.1.1 Precursor gas emissions

Sulfur emissions are dominated by SO${}_{\mathrm{2}}$ with 65452 $\pm$ 4483 $\mathrm{Gg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$ (model mean plus/minus standard deviation; see Table ). These emissions are largely of anthropogenic origin, as shown in Fig. , where we compare the natural-only versus all emissions in the NAT and REF simulations (from the MACCity inventory). Under anthropogenically perturbed conditions, emissions dominate in the Northern Hemisphere (NH) and are widely present over the ocean as well, in comparison to only natural emissions. In Fig. , ULAQ–CCM has slightly higher and WACCM6–CARMA lower SO${}_{\mathrm{2}}$ emissions than the other models. The second-strongest emissions are those of DMS, which show a large model spread, with a mean and standard deviation of 28169 $\pm$ 4453 $\mathrm{Gg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$. DMS emissions vary the most (with a coefficient of variation of 16 %; see Table ). Models typically calculate DMS emissions online as a function of wind speed and prescribed concentration of DMS in seawater. Large variations in DMS emissions have also been shown in a previous model comparison by . As an important precursor for tropospheric SO${}_{\mathrm{2}}$ (see fluxes provided by ECHAM5–HAM and SOCOL–AERv2 in Fig. ), this influences tropospheric SO${}_{\mathrm{2}}$ burdens. Figure  also shows the model results for OCS and H${}_{\mathrm{2}}$S for the models in which surface mixing ratios are specified instead of emissions. Their contribution to other sulfur-containing species can be tracked via the chemical fluxes provided by the SOCOL–AERv2 model. Their net surface flux is also shown for SOCOL–AERv2, calculated by balancing all other fluxes. Furthermore, ECHAM6–SALSA, WACCM6–CARMA, and WACCM6–MAM4 also include primary sulfate aerosol emissions. For example, in ECHAM6–SALSA, 2.5 $\mathit{\%}$ of all anthropogenic, wildfire, and volcanic sulfur emissions are emitted as sulfate aerosol particles. Finally, CS${}_{\mathrm{2}}$ is included by prescribing emissions; CS${}_{\mathrm{2}}$ is the dominant precursor to the important OCS (but not all models include this process in their chemistry scheme, since OCS is defined by a fixed mixing ratio). performed sensitivity experiments with the SOCOL–AERv1 model by isolating the contribution of specific sulfur emissions to the stratospheric aerosol layer. They found that OCS contributes around half to the stratospheric burden, SO${}_{\mathrm{2}}$ one-third, and the rest is mostly from DMS. However, their estimate for the total stratospheric aerosol burden was 109 $\mathrm{Gg}\left(\mathrm{S}\right)$, which is lower than the value obtained with the new model version (167 $\mathrm{Gg}\left(\mathrm{S}\right)$ in Fig. ).

3.1.2 Global tropospheric burdens

Although the implementation details of the sulfur cycle vary widely among the different models, most of them include the following three dominant species: OCS, SO${}_{\mathrm{2}}$, and sulfate aerosol. Calculating the model means of all species in Fig. , OCS accounts for about 68 $\mathit{\%}$ of the total mass of sulfur in the troposphere and varies very little among models with a coefficient of variation of 4 % (see Table ). This is not surprising because this species is prescribed at the surface and has few chemical sinks in the troposphere, while sinks such as plant uptake are not explicitly modeled. Sulfate and SO${}_{\mathrm{2}}$ represent 17 $\mathit{\%}$ and 10 $\mathit{\%}$ of the total tropospheric sulfur mass, respectively, though inter-model differences are much larger than for OCS. SO${}_{\mathrm{2}}$ burdens range from 108 $\mathrm{Gg}\left(\mathrm{S}\right)$ in CAM5–CARMA to 5 times that amount in ULAQ–CCM. Tropospheric aerosol also varies between the lowest values in ECHAM5–HAM at 312 Gg(S) and the highest in UM–UKCA (885 $\mathrm{Gg}\left(\mathrm{S}\right)$). However, in all models, these species are dominant and leave about 5 $\mathit{\%}$ tropospheric sulfur in the form of DMS, MSA, H${}_{\mathrm{2}}$S, CS${}_{\mathrm{2}}$, gas-phase H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, and SO${}_{\mathrm{3}}$. The latter species are present only in small amounts as a result of their respective short lifetimes (1 to a few days in the atmosphere) or low vapor pressure in the case of H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$. In contrast, SO${}_{\mathrm{2}}$ remains in the atmosphere for about 10 d and OCS for 2 years . However, the less abundant species fulfill an important role as precursor gases for sulfate particles. In addition, there are intermediate products of photolysis and oxidation and references therein, but they are even more short-lived and so are not shown here.

3.1.3 Global stratospheric burdens

The stratosphere shows a similar distribution of relative abundances as in the troposphere, with OCS accounting for about 64 % total stratospheric sulfur, while aerosol is 31 % and SO${}_{\mathrm{2}}$ is 2.5 % (see Table ). Based on its abundance and its lifetime, OCS was identified as the main precursor for background sulfate aerosol by (see also ). OCS is primarily a product of photolysis and subsequent oxidation of other precursor gases, such as CS${}_{\mathrm{2}}$ . It is characterized by a long tropospheric lifetime and is removed from the atmosphere mainly via plant uptake in the troposphere and photo-oxidation to H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ in the stratosphere . The variation between the models is 16 %, i.e., somewhat larger in the stratosphere than in the troposphere, with outliers at as much as 28 % above the model mean (ULAQ–CCM). The lifetime of OCS is still very long in the stratosphere (years). Therefore, the importance of dynamical (meridional and cross-tropopause transport) and chemical (photolysis) sinks increases. This increases the multi-model spread, given that individual models necessarily have some biases in their representation of dynamical processes. This is discussed further in the context of the spatial distribution in Sect. . Models tend to overestimate the stratospheric OCS burden when compared to MIPAS. The MIPAS dataset only extends to 35 km (Table ). However, not much OCS is present above this altitude due to photolysis (see Fig. ).

SO${}_{\mathrm{2}}$ represents an important step in the atmospheric sulfur oxidation chain, as all precursor gases are first oxidized to SO${}_{\mathrm{2}}$, from there to SO${}_{\mathrm{3}}$, and then to gaseous H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, which co-condenses with H${}_{\mathrm{2}}$O to form aerosol particles. The stratospheric SO${}_{\mathrm{2}}$ burden reveals the largest scatter among the three major species (coefficients of variation in Table  of $\pm$ 56 %), with a factor of 5 difference between the models with the smallest and largest burdens. Although the model spread is large, MIPAS is close to the overall mean of stratospheric SO${}_{\mathrm{2}}$ burden. We can group the models in Fig.  into three groups with respect to SO${}_{\mathrm{2}}$. SOCOL–AERv2, CAM5–CARMA, and UM–UKCA have the lowest burdens, while WACCM6–CARMA, WACCM6–MAM4, and MIROC–CHASER are in the middle range (and closest to MIPAS), and ULAQ–CCM, ECHAM6–SALSA, and ECHAM5–HAM have significantly higher burdens of SO${}_{\mathrm{2}}$ than the other models. The spatial distribution and plausible reasons for these differences are discussed in Sect. . Sources of SO${}_{\mathrm{2}}$ in the stratosphere include photo-oxidation of OCS and, just as important, the transport of tropospheric SO${}_{\mathrm{2}}$ across the tropopause. Uncertainties remain regarding the transport of SO${}_{\mathrm{2}}$ to the stratosphere , which, besides the complexity of the upper troposphere and lower stratosphere (UTLS) transport processes themselves, could also be very sensitive to how the models treat the aqueous oxidation processes in the upper troposphere . Tropospheric aqueous and gas-phase SO${}_{\mathrm{2}}$ oxidation fluxes are presented in Fig.  for SOCOL–AERv2 and ECHAM5–HAM; however, ECHAM5–HAM shows lower values for both pathways. In the lower stratosphere, SO${}_{\mathrm{2}}$ lifetime for oxidation by OH is 3–4 weeks, making it the source of H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ and subsequently the Junge aerosol. From these considerations, it appears that most of the models underestimate the turnover in the chemical reaction of SO${}_{\mathrm{2}}$ $+$ OH.

The stratospheric aerosol burden is about 156 $\pm$ 51 Gg sulfur in the multi-model mean (with standard deviation; see Table ). This value differs only slightly from the value derived from the SAGE-3$\mathit{\lambda }$ observational dataset described by and personal communication with Beiping Luo, 2023. The scatter between the models is rather small; however, it still has a coefficient of variation of 33 $\mathit{\%}$ with a factor of 3 difference between model outliers. The reduced scatter for aerosol compared to SO${}_{\mathrm{2}}$ is due to the much longer residence time of aerosol in the stratosphere, as well as the contribution of tropospheric aerosol. A slightly shorter or longer chemical lifetime of SO${}_{\mathrm{2}}$ does therefore not affect the SO${}_{\mathrm{2}}$ burden and the aerosol burden equally. The scatter of aerosol also does not match the scatter of OCS in the stratosphere because OCS accounts for only about half of the source of stratospheric aerosol, with the other half being the cross-tropopause fluxes of SO${}_{\mathrm{2}}$ and aerosol itself . In our analysis, unfortunately, only the net cross-tropopause fluxes are available and also only from a few models, which limits the analysis in the attribution of stratospheric biases. In addition, the model differences in stratospheric aerosol loading can also be caused by differences in sedimentation fluxes among models, which in turn depends on aerosol size distribution. and estimated that gravitational sedimentation reduces the background stratospheric aerosol burden by about half. In comparison, ECHAM5–HAM exhibits an unusually low aerosol burden. This is likely the result of anomalous concentrations of OH and, thus, slower SO${}_{\mathrm{2}}$ oxidation in the stratosphere, as discussed in Sect. . Other factors may also contribute, such as differences in the vertical residual velocity in the tropics, which was identified by as a cause of major differences in the aerosol burdens between ECHAM5–HAM and WACCM6. The highest aerosol values are reported by ECHAM6–SALSA. discussed how excessive new particle formation may contribute to such an effect. However, it is unclear if this nucleation bias only applies to scenarios with a large disturbance from volcanic eruptions or stratospheric aerosol injections. Conversely, this has not been observed before in ECHAM6–SALSA see, e.g.,, with the main difference being the vertical resolution of the model.

Short-lived species, such as DMS and gaseous H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ with lifetimes shorter than 1 d, exhibit large scatter between models (with models differing by more than 1 order of magnitude), but the uncertainties in these burdens do not significantly affect the ones of longer-lived species. Rather, one would need to investigate the reaction rates to determine the processes leading to differences in the burdens of the major species.

Figure  presents the total integrated sulfur burdens. The total emission and deposition rates in ECHAM6–SALSA do not include OCS but are in the same range as most other models, while the total S burden is lower by about the amount of OCS in the other models. However, ECHAM6–SALSA has the highest SO${}_{\mathrm{2}}$ burden out of all models, making more of it available for oxidation and the formation of aerosol. Deposition is again in the same range as other models, meaning that aerosol may accumulate more.

Figure 2

The total sulfur balance for each model. All sulfur species are summed up for each model and depicted here (in $\mathrm{Gg}\left(\mathrm{S}\right)$). Cross-tropopause fluxes (in $\mathrm{Gg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$) are net fluxes of aerosol precursor gases (net upward) plus the net aerosol flux (net downward), which were only provided by SOCOL–AERv2 and ULAQ–CCM. The burden for ECHAM6–SALSA is lower, as the model does not include OCS. Therefore, this model is excluded from the calculation of the model spread for the total sulfur burden.

[Figure omitted. See PDF]

3.1.4 Wet and dry deposition rates

Deposition rates are mainly available for sulfate aerosol, which is dominated by wet deposition, and for SO${}_{\mathrm{2}}$, which is dominated by dry deposition, consistent among all models. Some species, such as CS${}_{\mathrm{2}}$ or DMS, are often considered to be entirely chemically processed and therefore not deposited. While wet and dry deposition rates vary considerably (see for SO${}_{\mathrm{2}}$ or aerosol in Fig. ), the total deposition of sulfur remains similar among models, as seen in Fig. . This is not surprising because the total sulfur emission flux is meant to be the same (or at least similar) for all models. All emitted sulfur has to be deposited back to the surface, assuming that models do not have large errors in mass conservation. ECHAM5–HAM is a strong outlier in terms of total sulfur deposition, with the reason being currently unclear, given that other fluxes available for this model (DMS conversion to SO${}_{\mathrm{2}}$, SO${}_{\mathrm{2}}$ emission, and oxidation fluxes) do not differ much from other models. discussed the partitioning of sulfur deposition in a global model intercomparison for tropospheric aerosol. They found that, in some models, sulfur is already deposited in the form of precursor gases, therefore resulting in a lower burden and less deposition in the aerosol phase. Here, we see that ULAQ–CCM, ECHAM5–HAM, and WACCM6–MAM, which have lower tropospheric aerosol burdens, also exhibit lower total aerosol deposition rates compared to the other models. However, when adding all reported deposition fluxes together, ECHAM5–HAM has the highest percentage (79 $\mathit{\%}$ deposited sulfur) deposited as aerosol, while ULAQ–CCM has the lowest at 45 $\mathit{\%}$. Therefore, in ECHAM5–HAM another process must be influencing the deposition rates. Most aerosol is wet deposited, which is consistent with . The total aerosol deposition shown in Fig. a–f indicates that it is not just the amount but the spatial distribution of the aerosol deposition that varies considerably among models, though it still mostly resembles the global distribution of precipitation e.g.,. Regional differences are influenced by both the aerosol formation processes and the biases in the models' precipitation patterns and the details of deposition schemes . In all models, the anthropogenic aerosol (Fig. g–i) deposits mostly in the regions where the anthropogenic emission of SO${}_{\mathrm{2}}$ takes place (Fig. ). This confirms the short tropospheric lifetime of the aerosol, limiting its long-distance transport. In ECHAM5–HAM, 98 $\mathit{\%}$ of the aerosol deposition is wet deposition, followed by 95 $\mathit{\%}$ in ECHAM6–SALSA, 92 $\mathit{\%}$ in SOCOL–AERv2, 90 $\mathit{\%}$ MIROC–CHASER, 86 $\mathit{\%}$ in WACCM6–MAM4, and finally ULAQ–CCM with the lowest percentage at 82 $\mathit{\%}$ wet deposition.

For SO${}_{\mathrm{2}}$, the dry deposition dominates but varies strongly between models, with 99 $\mathit{\%}$ for SOCOL–AERv2, 92 $\mathit{\%}$ for ULAQ–CCM, 77 $\mathit{\%}$ for ECHAM5–HAM, 71 $\mathit{\%}$ for ECHAM6–SALSA, and only 59 $\mathit{\%}$ in MIROC–CHASER.

3.1.5 Total atmospheric sulfur burden

We present the total sulfur burden in the stratosphere and troposphere in Fig. . Although the models have chemistry schemes of various complexities, which affects the partitioning between the sulfur species, the total sulfur burden is similar across models, at 3861 $\pm$ 294 $\mathrm{Gg}\left(\mathrm{S}\right)$, for the models including the main three sulfur species OCS, SO${}_{\mathrm{2}}$, and sulfate aerosol (therefore excluding ECHAM6–SALSA). This corresponds to 485 $\pm$ 85 $\mathrm{Gg}\left(\mathrm{S}\right)$ in the stratosphere (12.6 $\mathit{\%}$ of all atmospheric sulfur) and 3375 $\pm$ 284 $\mathrm{Gg}\left(\mathrm{S}\right)$ in the troposphere. The relative difference among models is higher in the stratosphere, as discussed in Sect. . Figure  also presents estimates from SOCOL–AERv2 and ULAQ–CCM of the net cross-tropopause fluxes of aerosol and aerosol precursors. The two models agree that the sign of the net fluxes is directed upward for precursors and downward for aerosol but disagree on their magnitudes. These fluxes are not fully balanced, although they are expected to be somewhat in equilibrium. The reason for ULAQ–CCM is that only the SO${}_{\mathrm{2}}$ flux is available, and the contribution of other species (mainly OCS and DMS, as seen from SOCOL–AERv2 data in Fig. ) is missing. SOCOL–AERv2 takes all the species into account; therefore, the imbalance could be a result of deficient mass conservation in its transport scheme, especially where the tracer gradients are steep . As for emissions and deposition fluxes, the values provided in the figure are quite scattered because not all variables calculated by models were present in the output provided in the model database and could therefore not be included in our budget calculations. This also explains the inconsistency between the total deposition and the emission fluxes for individual models. All the required output was only provided for SOCOL–AERv2. As a result, the full emission and deposition fluxes for this model are calculated by including the chemical fluxes of minor precursors to SO${}_{\mathrm{2}}$, resulting in a good agreement between the two. In all other models, these data were not provided in the output, therefore statements about their mass conservation cannot be made. In terms of the model mean values presented in Table , individual model nuances are averaged out, resulting in a flux of 98 $\mathrm{Tg}\left(\mathrm{S}\right)$ for both the total emission and deposition.

3.2 Seasonal cycle of sulfur compounds

Figure  shows the seasonal cycle of the most abundant sulfur compounds in the atmosphere, as well as one of the main oxidizing agents of the atmosphere, the hydroxyl radical (OH). Also depicted in (Fig. a–e) is the relationship between the stratospheric sulfate burden and the stratospheric aerosol optical depth (SAOD) at $\sim$ 0.5 $\mathrm{\mu }$m for each model, season, and latitude. We see for all species, the multi-model spread is the smallest in the tropical region (here at 30° N–30° S), where the burdens are also lower. The relationship between SAOD and the aerosol burden shows little seasonal dependence and very little scatter. In the extratropics, the burdens of sulfur species are higher, suggesting additional transport through the subtropical tropopause and a larger spread among the models in terms of the relation of SAOD and aerosol mass. This suggests a larger divergence in the aerosol size distribution in the extratropics than in the tropics. Additionally, we see a distinct seasonal cycle in all species outside the tropical region and an increased seasonality of the SAOD/aerosol mass relation in many models.

Figure 3

The seasonal cycle of the three most abundant atmospheric sulfur species (sulfate aerosol, SO${}_{\mathrm{2}}$, and OCS), as well as OH (OH was only reported by five models). All data were zonally averaged over the latitude band given in each column title. Additionally, panels (a)–(e) represent the relationship between the SAOD at 0.5$\mathit{\mu }m$ and the sulfate burden for five latitude bins. Each season is depicted by a character, with a star for winter, a triangle for spring, a circle for summer, and a diamond for autumn. Panels (f)–(j) depict the seasonal cycle of the stratospheric aerosol burden for each of these latitude bins. In panels (k)–(o), we see the stratospheric OCS burden, while panels (p)–(t) show the stratospheric SO${}_{\mathrm{2}}$ burden and panel (u)–(y) the stratospheric OH burden, respectively. Each burden is given (in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{2}}$), except for OH (given as $\mathrm{kg}\left(\mathrm{OH}\right){\mathrm{m}}^{-\mathrm{2}}$). All data have been averaged over the 20 simulated years for each month. Not shown in the figure are the maximum values of the stratospheric SO${}_{\mathrm{2}}$ burden in ULAQ–CCM and the stratospheric OH burden in both WACCM6 models. ULAQ–CCM SO${}_{\mathrm{2}}$ burdens in the northern mid-latitudes peak in January at 1.9 $\times$ 10${}^{-\mathrm{7}}$ $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{2}}$ or 5.5 $\times$ 10${}^{-\mathrm{7}}$ $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{2}}$ in the northern polar region. Maximum OH burdens in WACCM6–CARMA and WACCM6–MAM4 in both the NH (July) and SH (December) amount to 1.3 $\times$ 10${}^{-\mathrm{8}}$ $\mathrm{kg}\left(\mathrm{OH}\right){\mathrm{m}}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

In addition to an agreement with Fig.  in terms of model aerosol burden levels, we observe that the models with higher burdens (ECHAM6–SALSA, ULAQ–CCM) have a stronger seasonal cycle than those with lower burdens (MIROC–CHASER, ECHAM5–HAM, and UM–UKCA). ECHAM6–SALSA has the highest burdens in the low and mid-latitudes. ULAQ–CCM is similar to ECHAM6–SALSA in the NH polar region and even exceeds the burdens in ECHAM6–SALSA in the Southern Hemisphere (SH). For aerosol, we observe a minimum at high latitudes in the vortex (Fig. f) from July to October, followed by an increase which can be attributed to the breakup of the southern polar vortex. In winter, the polar vortex edge represents a barrier to mixing and transport, preventing extra-polar sulfate from entering, similar to what has been observed for ozone and other species . In SAGE-3$\mathit{\lambda }$, this minimum occurs about 1 month earlier than in the models. The same pattern is seen for OCS in (Fig. k–o). We conclude that this effect is due to dynamics rather than chemistry. This minimum is discussed in more detail in Sect. .

In the northern mid-latitudes, most models and SAGE-3$\mathit{\lambda }$ have minimal aerosol burdens in June, July, and August. In contrast, ECHAM6–SALSA and ULAQ–CCM have increasing burdens during this time, as also seen in high northern latitudes. In the 60°–90° N latitude band, however, the seasonality is not as clear in the models as in the other latitude bands. ULAQ–CCM has two distinct peaks around May/April and again in September/October, while ECHAM6–SALSA and MIROC–CHASER have the highest burdens in this region only in October/November; ECHAM5–HAM and both WACCM6 models, as well as SAGE-3$\mathit{\lambda }$, can be grouped together with burdens peaking around May, while in SOCOL–AERv2, a seasonal cycle is not distinguishable. This indicates higher uncertainties in the models in this region, which is directly related to large differences in the northern polar vortex climatology and dynamics between the models , as also shown in Fig.  for the models participating here. In Fig. a.I–i.I in December–February (DJF), zonal winds are weaker and temperatures higher than ERA-I data in the northernmost latitudes in ULAQ–CCM, UM–UKCA, and to a lesser extent also in ECHAM6–SALSA above $\sim$40 hPa and ECHAM6-HAM and MIROC–CHASER below $\sim$40 hPa. These differences are, however, rather small in most models (with the largest outlier with temperature anomalies above 7 K being ULAQ–CCM), especially when comparing to differences in the SH, where seasonal cycles are very similar among models despite stronger biases in the southern polar vortex winds and temperatures, as seen in Fig. .

Another feature is the similarity between WACCM6–CARMA and WACCM6–MAM4. Considering that the two WACCM6 models share the same dynamical core and chemistry but are coupled to different aerosol microphysical models, this would indicate that the latter does not influence the resulting burden. There are some temperature differences among the two WACCM6 model configurations in the lower stratosphere; however, the statistical significance of this difference has not been tested. A comparison of the three ECHAM-based models (SOCOL–AERv2, ECHAM6–SALSA, and ECHAM5–HAM) is not as straight-forward, as they have different chemistry schemes and different vertical resolutions since they are based on different ECHAM versions, i.e., ECHAM5 in SOCOL–AERv2 and ECHAM5–HAM versus ECHAM6 in ECHAM6–SALSA. We also discuss, in Sect. , how an internally mixed type of aerosol (including other components additionally to sulfur) in CARMA causes differences in the troposphere. Furthermore, the effective radius varies, as discussed in Sect. , with potential impacts on SAOD.

The relationship between the SAOD and stratospheric aerosol burden is illustrated in Fig. a–e. (We use the pure sulfate SAOD for SOCOL–AERv2 and ULAQ–CCM, whereas in WACCM6–CARMA, it is the SAOD of both pure and mixed aerosol and in all other models, it is the total SAOD). In general, a higher stratospheric aerosol burden also leads to higher SAOD (Fig. a–e). This relationship is less clear in the extratropics. The relationship remains highly linear for SAGE-3$\mathit{\lambda }$, which is connected to the assumptions on size distributions in the construction of this dataset. As the relationship between burdens and SAOD is very similar among models in the tropics, the conversion schemes in all models are likely to be very similar. The differences in the extratropics hint at different size distributions, which influence the SAOD but not the total mass. This is discussed briefly in Sect. . Here, WACCM6–MAM4 tends to have lower SAOD than WACCM6–CARMA, revealing differences between the two aerosol schemes.

The seasonal variability in the OCS burden is very similar to that of sulfate aerosol (Fig. k–o). The spread in burdens and irregularity in the seasonal cycles appears slightly smaller for OCS. The highest scatter is found in the northern high latitudes, with CAM5–CARMA now also displaying a second peak, as seen in ULAQ–CCM, around October. In MIPAS, the minimum in the southern polar region is about a month earlier than in all models.

For SO${}_{\mathrm{2}}$ in Fig. p–t, we can again group models with low and high burdens, as specified in Sect. , with ULAQ–CCM, ECHAM6–SALSA, and ECHAM5–HAM in the higher range; the WACCM6 models and MIROC–CHASER in the mid-range; and SOCOL–AERv2, CAM5–CARMA, and UM–UKCA at the lower end. Possible reasons for higher SO${}_{\mathrm{2}}$ burden could be less oxidation or more transport across the tropopause. The stratospheric OH burden is provided in Fig. u–y. In ECHAM5–HAM, the OH burden is distinctly lower, supporting the hypothesis that less SO${}_{\mathrm{2}}$ is oxidized and hence less H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ is available for aerosol formation. have previously shown how the prescribed OH fields in ECHAM5–HAM change the aerosol burden after the 1815 Tambora eruption compared to models with interactive OH chemistry. In ULAQ–CCM, stratospheric OH is very close to the values of other models. Additionally, while ECHAM5–HAM is among the highest burdens in all latitudes, burdens in ULAQ–CCM are only elevated in the polar regions and northern mid-latitudes. Figure  shows the temperatures in each model compared to the 1999-2004 ERA-I average. In ULAQ–CCM in Fig. b.I–b.IV, the winter polar stratosphere is warmer and has weaker westerlies compared to ERA-I, indicating a weaker polar vortex and hence a weaker barrier to transport from mid-latitudes to the poles. The distribution of sulfate and SO${}_{\mathrm{2}}$ in ULAQ–CCM will be discussed in more detail in Sect. .

Stratospheric OH is too short-lived to be transported. It is mainly formed by the reaction of H${}_{\mathrm{2}}$O with O(${}^{\mathrm{1}}$D), which is a product of O${}_{\mathrm{3}}$ photolysis . Therefore, less OH is available in the wintertime in the middle and high latitudes. Since OH is the most important oxidizing agent for SO${}_{\mathrm{2}}$, we expect lower OH burdens to correspond to higher burdens of SO${}_{\mathrm{2}}$. We see that for those models where OH was provided, SOCOL–AERv2, ULAQ–CCM, and WACCM6–CARMA have higher stratospheric OH loadings than ECHAM5–HAM, with the latter also having a much less pronounced seasonal cycle. However, these differences do not directly translate to differences in SO${}_{\mathrm{2}}$ burden, as shown by ULAQ–CCM, which has a higher SO${}_{\mathrm{2}}$ burden than the other models. discussed the importance of interactive OH during volcanically perturbed conditions, where OH is depleted more rapidly by SO${}_{\mathrm{2}}$ oxidation. In non-interactive chemistry schemes, such as in ECHAM5–HAM and ECHAM6–SALSA, the OH fields may need to be adapted for the conditions, or there could be too much OH available after volcanic eruptions if the background OH level is used but too little during quiet conditions if the volcanically depleted OH is used.

3.3 Spatial and vertical distribution

3.3.1 Sulfate aerosol

Figure  shows the distribution of aerosol in all models (Fig. a–i) and SAGE-3$\mathit{\lambda }$ (Fig. k). The relative standard deviation (RSD, expressed as a percentage of the multi-model mean) distribution in Fig. l indicates the regions where the inter-model differences are the largest. Most of the aerosol mass is concentrated close to the surface, especially in tropical regions, and more pronounced in low northern latitudes. This is likely due to higher SO${}_{\mathrm{2}}$ concentrations in these regions, where anthropogenic emissions originating from East Asia are dominant . Figure g–l compares the sulfate mass mixing ratios of the NAT experiment with the REF experiment. It shows how anthropogenic emissions increase aerosol concentrations by 30 % and up to more than 80 % locally in the NH troposphere and lower stratosphere. The SH is less influenced by anthropogenic emissions. The same feature is seen in the sulfate deposition in Fig. g–l, which increases by more than 80 % over large parts of the NH continents. In SOCOL–AERv2 and ULAQ–CCM, the influence of anthropogenic emissions extends higher into the stratosphere, while in ECHAM5–HAM and MIROC–CHASER, this additional SO${}_{\mathrm{2}}$ and aerosol is likely removed before reaching the stratosphere. The WACCM6 models retain some of their similarities discussed in the previous section in Fig. . Differences can be seen in the lower troposphere, where mixed aerosols (WACCM6–CARMA) gain importance. However, this does not apply to CAM5–CARMA, which resembles WACCM6–MAM4 more in the troposphere, despite sharing the microphysical scheme with WACCM6–CARMA. The RSD increases rapidly above the lowest model layers in the tropical region, highlighting possible model differences in tropical upwelling and removal of aerosol by convective precipitation. Higher values of the RSD of about 50 % also extend latitudinally in the subtropical tropopause region, indicating model differences in the UTLS transport. Very high values are also seen in the extratropical troposphere. The best model agreement is found in the Junge layer region. The high values above 10 hPa can be disregarded, as very little aerosol is expected to be at this height.

Figure 4

Stratospheric aerosol mass density (in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{3}}$) and averaged over the 20 simulation years (10 years for ULAQ–CCM) for each of the models (a–i) in comparison with 6 quiet years from 1999–2004 for SAGE-3$\mathit{\lambda }$ (k). The dark blue lines are the time-averaged tropopauses in each respective model (a–i) and the ERA-I tropopause, using the World Meteorological Organization definition for SAGE-3$\mathit{\lambda }$ (k). To obtain the mass density, we converted the mass mixing ratios of sulfate aerosol using the ideal gas law and provided temperature and pressure fields. Panel (l) shows the relative standard deviation as a percentage of the multi-model mean, where all model data were interpolated to 39 pressure levels and gridded onto a 5° $\times$ 5° grid.

[Figure omitted. See PDF]

In SAGE-3$\mathit{\lambda }$, the elevated values in the tropical troposphere are likely caused by minor volcanic eruptions. mention four volcanic eruptions during 1999 and 2004: Ulawun (September 2000), Shiveluch (May 2001), Ruang (September 2002), and Reventador (November 2002). This last eruption emitted the largest amount of SO${}_{\mathrm{2}}$, with estimates by at around 84 Gg, and a plume height of about 17 km or 94 Gg, reaching a plume height of about 22 km .

Figure 5

The stratospheric aerosol column (in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$), averaged for each month over the 20 years of simulation, 10 years for ULAQ–CCM, and 6 years for SAGE-3$\mathit{\lambda }$. (a–i) Each panel represents one model for all latitudes over time. Panel (k) is the SAGE-3$\mathit{\lambda }$ dataset. Panel (l) shows the RSD of all models (not including SAGE-3$\mathit{\lambda }$). All data were gridded to a 5° $\times$ 5° grid to calculate the RSD.

[Figure omitted. See PDF]

Considering only the stratosphere in Fig. , maxima are seen in the mid-latitudes and polar regions around the tropopause. Since emissions in the NH exceed those in the SH , and the stratospheric meridional transport from the tropics is weaker in the SH, the overall SH aerosol burden is also slightly lower in most models and in SAGE-3$\mathit{\lambda }$. Dynamical processes in the UTLS, such as isentropic mixing may have a strong influence on the stratospheric burden. The strongest maximum is seen in ULAQ–CCM (see also Fig. ), with most of the aerosol residing above the tropopause. In all models, in the northern subtropical to mid-latitudes, the higher tropospheric burdens extend to higher altitudes, even “connecting” to the stratosphere. argued that 15 $\mathit{\%}$ of NH stratospheric aerosol originates from the Asian summer monsoon, where anthropogenic emissions are high (see also Fig. ). Additionally, this could also be a result of higher emissions in the NH, coupled with generally higher convection above the continents and isentropic transport and mixing with stratospheric air . As ULAQ–CCM has the lowest horizontal resolution, which leads to weakened transport barriers and stronger numerical diffusion, this effect could be a defining factor in this model. Another factor that could be responsible for the model biases is the climatology and variability in the tropopause. Figure  shows a high multi-model spread of the tropopause pressure, especially at high latitudes.

The seasonal cycle of the stratospheric aerosol burden is shown in Fig. . The main outliers in terms of absolute values are ULAQ–CCM, ECHAM6–SALSA, and ECHAM5–HAM, which is consistent with the values in Fig. . As in Fig. , the tropical region is marked by low values without an obvious seasonality. This is also evident from the RSD in Fig. l, with the lowest values in the tropics. Mid-latitude burdens coincide with the seasonal tropopause shift, where there is a sharp gradient towards the pole around 30°–50° N and S. In the SH, the highest burdens extend the furthest north to about 40° S in austral winter around August and September. In the NH, the pattern is not as smooth, and the RSD is also generally larger than in the SH. The largest RSD values lie in the northern mid-latitudes, while most other areas have an RSD of about 20 $\mathit{\%}$–40 $\mathit{\%}$. As already described in Fig. , this mirrors the higher tropospheric burdens at this latitude in ULAQ–CCM and is situated at the same time and place as the Asian summer monsoon . All models and SAGE-3$\mathit{\lambda }$ show an aerosol minimum within the southern winter polar vortex due to the transport barrier in the winter and dominant downward transport within the vortex that brings sulfur-poor air from above. The RSD is higher in this region, possibly due to inter-model differences in vortex isolation.

Figure 6

The OCS mass density (in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{3}}$) is shown in the top row (a–c). The bottom row depicts the zonally averaged stratospheric OCS burden (in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{2}}$) (d–f). Panels (a) and (d) show the MIPAS dataset, while panels (b) and (e) are the multi-model means, and panels (c) and (f) are the RSD. All model data were interpolated onto the same grid and vertical coordinates. The data were calculated from volume mixing ratios of OCS using the ideal gas law and provided temperature and pressure fields. ECHAM6–SALSA does not track OCS and is therefore not included here.

[Figure omitted. See PDF]

As the latitudinal distribution of OCS in Fig.  also shows a minimum within the SH polar vortex, we conclude that the model scatter is due to differences in the dynamics among models. The magnitude of the aerosol minimum varies little among models (Fig. l) and is also present in SAGE-3$\mathit{\lambda }$. Still, the timing of the vortex formation and breakup varies among models, and this has repercussions for aerosol transport. showed that the onset of the stratospheric final warming (SFW) tends to take place too late in many CMIP6 models in both hemispheres. WACCM- and MIROC-based models were part of this study, where the former tended to have a delayed SFW by up to 20 d, while MIROC tended to be much closer to the Japanese 55-year Reanalysis (JRA-55; ). In Fig. , we see that in WACCM6–CARMA and WACCM6–MAM4, SOCOL–AERv2, and CAM5–CARMA, for example, the aerosol is transported into the Antarctic region only in December, while this already takes place around November in ECHAM5–HAM, ULAQ–CCM, and UM–UKCA. In MIROC–CHASER, the transport seems to increase gradually over a longer period compared to other models, with values already increasing around September–October. In SAGE-3$\mathit{\lambda }$, the polar vortex aerosol minimum appears seasonally earlier in comparison to all models. From Fig. , which compares the zonal winds and temperature fields from the models to ERA-I data, it is evident that the southern winter polar vortex in WACCM6–CARMA and WACCM6–MAM4 is too strong in September–November (SON), while it is too weak (and warm) in ULAQ–CCM and ECHAM5–HAM. For the latter, this indicates an early onset of the SFW as opposed to WACCM6, where it takes place at a later time.

Larger differences are seen in the northern high to mid-latitudes, as already described in Sect. . Similar to the southern polar vortex, though much less pronounced, the northern polar vortex is marked by a local minimum in some models. This is only clear in ECHAM6–SALSA, CAM5–CARMA, WACCM6–CARMA, and WACCM6–MAM4. Although most models do depict a weak minimum, it is not well resolved in Fig. . From Fig. , we conclude that the northern polar vortex is too weak in ULAQ–CCM, while the signal is not as obvious in other models. A slightly stronger vortex is, for example, seen in ECHAM6–SALSA in December–February (DJF), while CAM5–CARMA, ECHAM5–HAM, and both WACCM6 models are slightly too cold in the lower and middle stratosphere. From March to May (MAM), the latter two have more pronounced northern polar vortices, indicating a generally delayed SFW, as with the southern polar vortex. For aerosol, this means too little transport into the polar region, resulting in a lower polar and higher mid-latitude optical depth.

3.3.2 Carbonyl sulfide

The vertical distribution of OCS is very similar in all models and the MIPAS observational data, with high burdens close to the surface and a relatively uniform gradient with height and much weaker latitudinal variation than for sulfate aerosol, as seen in the multi-model mean in Fig. . Averaging kernels were not considered here, which can improve the agreement of observations and model data, as shown in . In the tropics, OCS is found at slightly lower pressures than at high and mid-latitudes, due to the vertical transport in this region . Subsequently, OCS is oxidized and forms SO${}_{\mathrm{2}}$ at around 20 km (or ca. 60 hPa). Oxidation via O and OH only plays a minor role . The low inter-model spread is also reflected in the RSD, with very low values in most of the troposphere. As the mass of OCS decreases, the RSD increases with height. In Fig. b, we see that the largest RSD values are in the stratosphere; although, at altitudes above 10 hPa, this can be disregarded as not much OCS is present there.

In terms of seasonality of the stratospheric burden, OCS showed similar issues to the ones that were discussed for sulfur in individual models in the previous section, though the RSD is generally smaller. It is the largest in the polar regions, which is strongly related to the differences in the tropopause (see Fig. ). The tropopause in ULAQ–CCM extends lower in the polar regions than in other models, which with this uniform distribution of OCS allows more of it to reside in the stratosphere. This results in a pattern similar to sulfate in Fig. , where ULAQ–CCM has high burdens in the lower extratropical stratosphere, explaining its high total OCS burden in Fig. . Compared to MIPAS, the models behave very similarly, whereas again the minimum in the SH polar vortex is earlier and more pronounced in the observations compared to the multi-model mean.

3.3.3 Sulfur dioxide

SO${}_{\mathrm{2}}$ mass density is largest in the lower troposphere and the UTLS and decreases in the higher stratosphere (Fig. ). In some models, higher SO${}_{\mathrm{2}}$ values can also be seen in the tropical free troposphere. SO${}_{\mathrm{2}}$ burdens are higher directly below rather than above the tropical tropopause, which is opposite to the behavior of the aerosol burden in Fig. , as a result of the conversion of SO${}_{\mathrm{2}}$ to H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$, and subsequently sulfate, during the vertical transport of air in this region. Similar to Fig. , higher SO${}_{\mathrm{2}}$ values reach up to the tropopause in most models north of about 30° N. This is again a result of anthropogenic SO${}_{\mathrm{2}}$ emissions in this area, as shown in Fig. a–f. In REF, more SO${}_{\mathrm{2}}$ is transported to the stratosphere in many models. In ULAQ–CCM, WACCM6–CARMA, and WACCM6–MAM4, for example, the lowermost stratospheric SO${}_{\mathrm{2}}$ mixing ratio is increased by 50 % to up to 80 %, due to anthropogenic emissions. In contrast, in ECHAM5–HAM and MIROC–CHASER, most of these emissions do not reach the stratosphere. This then results in the inter-hemispheric asymmetry between the models in terms of the lower stratospheric SO${}_{\mathrm{2}}$.

Figure 7

The stratospheric SO${}_{\mathrm{2}}$ burden (in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{3}}$) on a logarithmic color scale. The models are depicted (a–i), while the MIPAS-derived dataset is shown (k). The RSD (l) is obtained from the multi-model ensemble. The tropopause is shown as a blue line and was averaged over the simulated time period and all longitudes. The mass density was obtained from volume mixing ratios, using the model temperature and pressure applied to the ideal gas law.

[Figure omitted. See PDF]

Figure 8

The stratospheric SO${}_{\mathrm{2}}$ burden (in $\mathrm{kg}\left(\mathrm{S}\right){\mathrm{m}}^{-\mathrm{2}}$) for all models (a–i), MIPAS (k), and the multi-model RSD (l). See also Fig. .

[Figure omitted. See PDF]

Although MIPAS (Fig. 7k) has a distribution more similar to, e.g., ECHAM6–SALSA than the models with lower stratospheric SO${}_{\mathrm{2}}$ burdens, the values in Fig.  are closer to, e.g., WACCM6–CARMA. In the stratosphere, most models and MIPAS have a secondary maximum of SO${}_{\mathrm{2}}$ in the tropical region between about 10-20 hPa. This can be attributed to the photolysis of OCS, discussed in the previous section, as well as to the evaporation of aerosol and subsequent photolysis of H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ . In MIPAS, we see an increase in SO${}_{\mathrm{2}}$ above 40 hPa, even outside the tropics, which is not seen in any of the models to the same extent. This is attributed to H${}_{\mathrm{2}}$SO${}_{\mathrm{4}}$ photolysis, which could be underestimated in models . As ECHAM6–SALSA does not include OCS, SO${}_{\mathrm{2}}$ only resides in the lower stratosphere. The other extreme is seen in ECHAM5–HAM, where too little SO${}_{\mathrm{2}}$ may be oxidized due to low OH concentrations (see Sect. ). Therefore, the lifetime of SO${}_{\mathrm{2}}$ increases, and there is more of it available to be further transported toward the poles from the tropics. The RSD in Fig. 7l is the highest in the UTLS and the northern extratropical troposphere, indicating low agreement among models in these regions. Models agree the most in the free tropical troposphere and the lower tropical stratosphere. Above 30 hPa, the RSD again increases, indicating differences in how models behave at the upper edge of the Junge layer both in terms of microphysics and photolysis.

Figure  shows the stratospheric SO${}_{\mathrm{2}}$ burden distributions. As for OCS in Figs.  and , the models follow a similar seasonal pattern in which the tropical region can be clearly distinguished from the extratropics. However, for SO${}_{\mathrm{2}}$, the lowest burdens are not found in the tropical region for all models. Instead, in CAM5–CARMA, the SH has lower burdens, while in most other models, these lower burdens are mainly seen in polar summer (except for ULAQ–CCM and UM–UKCA). Lower burdens are seen in the SH for most models, though still often exceeding the tropical burdens for several months, as in ECHAM5–HAM, WACCM6–MAM4, WACCM6–CARMA, and MIROC–CHASER. As described above, the higher burdens in ECHAM5–HAM, ECHAM6–SALSA, and ULAQ–CCM are found in the extratropical region. In WACCM6–CARMA and WACCM6–MAM, whose stratospheric SO${}_{\mathrm{2}}$ burdens were in the middle range, only the NH is characterized by higher burdens. While the higher SO${}_{\mathrm{2}}$ burdens in ECHAM5–HAM were previously explained by lower OH burdens, the higher SO${}_{\mathrm{2}}$ burdens in ECHAM6–SALSA and ULAQ–CCM can now be explained by increased transport across the tropopause in the northern subtropical region. While the differences among models are considerable, the best agreement outside the tropics is found from November to December at southern high latitudes (Fig. l).

Figure 9

The effective radius, SAD, and extinction for each model in the tropical region (30° N to 30° S) and in a single grid cell at Laramie (41° N, 105° W). Additionally, we show the effective radius and surface area density for SAGE II (observations in panels (a) and (b) and the extinction from GloSSACV2.2 in panel (c)). Observations in panels (d) and (e) are derived from OPC measurements. Models with sectional aerosol schemes are shown as solid lines, while modal aerosol schemes are dash-dotted lines. Extinction in SOCOL–AERv2, ULAQ–CCM, and MIROC–CHASER is at 525 nm; in CAM5–CARMA at 532 nm; in ECHAM6–SALSA at 533 nm; and in ECHAM5–HAM, WACCM6–CARMA, WACCM6–MAM4, and UM–UKCA at 550 nm. We use the effective radius of pure sulfate for SOCOL–AERv2 and ULAQ–CCM, an SAD-weighted average of pure sulfate and mixed aerosol effective radius for CAM5–CARMA and WACCM6–CARMA and the effective radius of all aerosol for ECHAM6–SALSA, ECHAM5–HAM, WACCM6–MAM4, and MIROC–CHASER. SAD describes only pure sulfate in SOCOL–AERv2, ULAQ–CCM, and CAM5–CARMA, while it encompasses all aerosol in ECHAM6–SALSA, ECHAM5–HAM, WACCM6–MAM4, and MIROC–CHASER. In WACCM6–CARMA, SAD is summed up for pure sulfate and mixed aerosol. Extinction of pure sulfate is used in SOCOL–AERv2 and ULAQ, while it is summed up for pure sulfate and mixed aerosol in WACCM6–CARMA, and for all remaining models, it is the total extinction.

[Figure omitted. See PDF]

3.4 SAD and effective radius

Figure a and d show the simulated effective radius, surface area density (SAD; Fig. b and e), and extinction in the tropics in (Fig. c), averaged over the whole time period; for the first two aerosol variables, values measured at Laramie, Wyoming (41° N, 105° W), are also shown. In the models, the closest grid box to Laramie was selected for comparison to the balloon observations. The effective radius in Fig. a increases above 100 hPa and reaches a maximum at around 20–30 hPa, ranging from about 0.12 $\mathrm{\mu }\mathrm{m}$ (UM–UKCA) to 0.22 $\mathrm{\mu }\mathrm{m}$ (ULAQ–CCM). At Laramie, this maximum is found at slightly lower altitudes at 40–60 hPa in all models, except ECHAM6–SALSA and ULAQ–CCM. While ULAQ–CCM and UM–UKCA are still outliers, all other models are closer together in this region in terms of their effective radius. In Sect. , we suggest that the size distribution is the cause of differences in SAOD between WACCM6–CARMA and WACCM6–MAM4, despite very similar burdens. Figure  confirms that the effective radius is larger in WACCM6–MAM4 (below ca. 15 hPa), while its extinction is lower than in WACCM6–CARMA. This is also seen at Laramie in Fig. d, although to a lesser extent. The SAGE-II-derived data indicate a smaller effective radius than in the models, except for MIROC–CHASER at this altitude, while exceeding the models' effective radius at lower levels. However, we have to stress that did not use the SAGE II data below 21 km (or ca. 50 hPa), due to the lower quality of the data at lower altitudes. Additionally, when the aerosol is low, particularly in background conditions, the SAGE-II-derived data and optical particle counters deviate significantly from each other . When comparing to the OPC measurements at Laramie, we also see a much better agreement with the models than with the SAGE II data in the tropical region. showed a distribution of effective radii in the pre-Pinatubo quiet conditions ranging from 0.1 $\mathrm{\mu }\mathrm{m}$ in SOCOL–AERv2 to 0.17 $\mathrm{\mu }\mathrm{m}$ in ECHAM6–SALSA and 0.27 $\mathrm{\mu }\mathrm{m}$ in ULAQ–CCM. In our study, however, ECHAM6–SALSA has the smallest effective radius out of these three models. In direct comparison with the volcanically perturbed values in , we see that the maximum effective radius is reduced to about half during these quiet conditions.

SAD is marked by a larger disagreement between the models, reflecting the differences among burdens described in Sect. . At 70 hPa, SOCOL–AERv2, ULAQ–CCM, WACCM6–MAM4, and MIROC–CHASER are closest to the observation-derived data by SAGE II, although the latter shows a sharper decrease towards 60 hPa, where WACCM6–CARMA, WACCM6–MAM4, and CAM5–CARMA are closer. Again, more models are closer to the OPC observations in mid-latitudes in Fig. 9e than to observations in Fig. 9b. In Fig. 9b, SAD in CAM5–CARMA and WACCM6–MAM4 decreases rapidly around 80 hPa, consistent with lower aerosol burdens in the upper tropical troposphere. In CAM5–CARMA, the SAD decreases to 0, as this variable is not defined below the tropopause.

Differences among models are again very large, concerning both the altitude and magnitude of the maximum extinction in the stratosphere (see Fig. c). While in CAM5–CARMA the maximum extinction is slightly below the tropopause (at ca. 100 hPa), it is above the tropopause for most other models. In MIROC–CHASER again, the stratospheric values are very low throughout and only increase in the troposphere. The largest stratospheric extinction coefficients are seen in ECHAM6–SALSA, which is consistent with the higher burdens in this model, followed by SOCOL–AERv2 and ULAQ–CCM. CAM5–CARMA and both WACCM6 models have mid-range values, consistent with GloSSACV2.2, while ECHAM5–HAM has the lowest extinctions, also consistent with its overall aerosol burdens. In , the extinction coefficients also vary a lot among models for volcanically perturbed conditions.

4 Conclusions

In this study, we present an analysis of the atmospheric sulfur budget using nine state-of-the-art global circulation models with interactive aerosol modules and sulfur chemistry. In particular, we showcase the burdens of the main sulfur species represented in these models (OCS, sulfate aerosols, and SO${}_{\mathrm{2}}$), as well as their latitudinal and vertical distribution and seasonal variability. Emission and deposition fluxes are also discussed. Sulfate aerosol, OCS, and SO${}_{\mathrm{2}}$ behave very differently in the atmosphere, providing insight into the modeled processes. OCS is almost chemically inert in the troposphere and is primarily influenced by transport mechanisms but undergoes photolysis in the stratosphere. Conversely, SO${}_{\mathrm{2}}$, which has a lifetime of a few days to a few weeks, is more responsive to chemical processes in the troposphere.

Previous studies have revealed the main sulfur species abundant in the atmosphere to be OCS, SO${}_{\mathrm{2}}$, and sulfate aerosol, which we here find to make up about 95 $\mathit{\%}$ of the total tropospheric and 98 $\mathit{\%}$ of the stratospheric sulfur mass, averaging over all models which include these three species. However, variability across models in the stratospheric burden of each of these main species is large, with burdens of 319 $\pm$ 50 $\mathrm{Gg}\left(\mathrm{S}\right)$ for OCS, 12 $\pm$ 7 $\mathrm{Gg}\left(\mathrm{S}\right)$ for SO${}_{\mathrm{2}}$, and 156 $\pm$ 51 $\mathrm{Gg}\left(\mathrm{S}\right)$ for sulfate aerosol. Anthropogenic emissions of SO${}_{\mathrm{2}}$ dominate the sulfate aerosol burden in the NH troposphere and increase the aerosol mass in the lowermost stratosphere by 10 %–80 %, which is very model-dependent. The SH, both in the troposphere and the stratosphere, is less affected by anthropogenic activity (10 %–30 %). The total deposition varies across models by about a factor of 2 between 99 692 and 51 211 $\mathrm{Gg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$ (or 86 693 $\pm$ 18 115 $\mathrm{Gg}\left(\mathrm{S}\right){\mathrm{yr}}^{-\mathrm{1}}$). Models also disagree on the form in which sulfur is deposited, as well as the split between wet and dry deposition. More sulfur being deposited as aerosol precursor gases could lead to a lower aerosol burden, as discussed by . We do, however, not explicitly observe such behavior. Furthermore, the relative importance of wet and dry deposition is expected to depend on the solubility of the compound, as well as the treatment of clouds and precipitation in the respective models . Evaluating model performances with respect to clouds and precipitation is beyond the scope of this study, and we refer to each model's documentation for details, as well as to dedicated model intercomparison activities on this topic e.g.,.

In terms of the spatial distribution, the models agree very well on the distribution and seasonality in the SH for all species. However, in the NH, where anthropogenic emissions contribute significantly to the aerosol burden, especially at low to mid-latitudes, more uncertainties persist. This is in part due to dynamical uncertainties in the northern polar vortex affecting the transport into high latitudes. The highest values in the northern extratropics are seen in ULAQ–CCM, which are associated with a slower and warmer vortex compared to the other models. Additionally, resolution-dependent isentropic transport through the subtropical tropopause could play a major role in increasing the extratropical burdens of OCS and sulfate in some models, which may be even more pronounced during the Asian monsoon, during which aerosol transport to the stratosphere is facilitated . A low horizontal resolution also increases diffusion across transport barriers, such as the polar vortex or the tropopause, where it crosses the isentropes in the subtropics . The speed of the BDC also could be a major factor of uncertainty affecting both the transport through the shallow branch, as well as upwelling, and thus the confinement of aerosol in the tropics. Idealized tracer experiments using SO${}_{\mathrm{2}}$ emissions might help to isolate the dynamical biases.

The analysis on the effective radius, SAD, and extinction in the tropics and at Laramie reveals notable differences among models. While the effective radius shows some agreement among models and is within the uncertainty in the OPC measurements in Laramie in all models, SAD exhibits larger discrepancies, reflecting the variations in aerosol burdens. Extinction levels also vary significantly, with ECHAM6–SALSA showing the highest values. Elaborating further on the reasons for the size distribution uncertainties would require a closer look at the individual microphysical processes, as was also highlighted by , but such data were not available in most of the models in our study. Recently, concerns have arisen regarding nucleation schemes, such as , which may overestimate the nucleation of new particles . However, these studies focus on stratospheric aerosol injection scenarios, where sulfur is much more abundant than in the background conditions presented in this study. Using observational data, showed that small volcanic eruptions affect the background size distribution of aerosols in unexpected ways, sometimes decreasing instead of increasing the effective radius. The volcanic events highlighted there can be further used as model test cases, contributing knowledge to the model performance in terms of microphysics and its dependence on the background conditions.

With this study, we have provided a first multi-model analysis of the background burdens of sulfur in atmospheric models. We mainly assessed the burdens and latitudinal and vertical distributions of sulfur species, allowing us to identify potential reasons for model diversity. As burdens are merely the result of fluxes, the next step would be to compare chemical and microphysical fluxes, as well as emission, cross-tropopause, and deposition fluxes. For SO${}_{\mathrm{2}}$ in particular, cross-tropopause fluxes in the mid-latitudes have the potential to provide more insight into the underlying processes causing differences in stratospheric burdens.

Overall, this study sheds new light on the atmospheric sulfur budget and, for the first time, presents it in a multi-model context. Our results highlight the importance of an interplay between the chemical, microphysical, and dynamical processes in atmospheric models. The interdependence of these processes complicates the attribution of biases when looking only at their final products, namely the distribution and variability in sulfur species in the atmosphere. The same holds true for the recent similar studies focused on volcanic events e.g., and stratospheric aerosol injections e.g.,. Some conclusions could be drawn by relating different variables. However, more detailed investigations of the reasons for inter-model disagreement would require process-oriented experiments such as the sensitivity experiments proposed by ISA–MIP for the Pinatubo eruption in the Pinatubo Emulation in Multiple models (PoEMS) e.g.,. However, these experiments still include the full range of complexity and may therefore be difficult to interpret. Much more illuminating would be experiments in which the degrees of freedom are reduced so that specific questions can be addressed. Comparison of standalone parts of the respective model codes, i.e., box model versions of the microphysical and chemical schemes, short and long wavelength radiative transfer modules, and treatment of the Mie calculation, is one possible solution. The other would be to perform dedicated tracer studies or simple sulfur cycle experiments (with a prescribed aerosol size distribution). Revisiting the ISA–MIP experiments in this sense and developing a set of dedicated benchmark simulations related to available observations would therefore be a valuable next step to validate and improve global stratospheric aerosol models but also to reduce uncertainties in solar radiation modification scenarios.

Appendix A

Figure A1

The SO${}_{\mathrm{2}}$ emissions in the NAT (a) and REF (b) experiments. NAT is composed of volcanic outgassing (c) and biomass burning (e), while the REF experiment additionally includes anthropogenic emissions (d).

[Figure omitted. See PDF]

Figure A2

Total deposition fluxes of sulfate aerosol (as a sum of dry and wet deposition) as absolute values are given in panels (a)–(f) and as a percent difference between the REF and NAT experiments in panels (g)–(l). The percent difference is calculated as REF–NAT and divided by REF. The data are time-averaged over the whole time period.

[Figure omitted. See PDF]

Figure A3

The difference in the time slice simulations averaged over 20 years for each season for temperature and zonal winds. Anomalies are calculated with respect to 1999–2004 ERA-I data, which was regridded to the respective model grid. Temperature is shown in color, while winds are in gray contours. Solid lines indicate stronger-than-observed westerly zonal winds, while dashed lines are weaker than observed. The numbers along with the contours are the anomalies (in $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$). Each column represents one season, while each row represents one model. The simulated tropopause is shown as a black line, also averaged over each season, respectively.

[Figure omitted. See PDF]

Figure A4

The percent difference between the REF experiment and the NAT experiment of the SO${}_{\mathrm{2}}$ volume mixing ratio in panels (a)–(f) and the sulfate mass mixing ratio in panels (g)–(l). The percent difference is calculated as REF–NAT divided by REF. All data are zonally averaged and over the whole time period of each simulation.

[Figure omitted. See PDF]

Figure A5

The tropopauses of each model, as well as the ERA-I-derived tropopause, which was used to calculate stratospheric burdens for SAGE-3$\mathit{\lambda }$.

[Figure omitted. See PDF]

Data availability

All model data used in this study can be accessed through 10.5281/zenodo.10411013 . The MIPAS combined SO${}_{\mathrm{2}}$ dataset is publicly available (10.35097/1788, ), as is MIPAS OCS (10.35097/1821, ). The GloSSACv2.2 data are available at 10.5067/GLOSSAC-L3-V2.2 . The SAGE-3$\mathit{\lambda }$ data are also based on this dataset. The ERA-Interim reanalysis is available from the Copernicus Climate Change Service (2023) (10.24381/cds.f2f5241d). OPC balloon measurements can be accessed through 10.15786/c.6379371. SAGE II v7 data are available from 10.5067/ERBS/SAGEII/SOLAR_BINARY_L2-V7.0 . The full raw dataset for all models is also available on the DKRZ Levante HPC system and can be made available upon request.

Author contributions

The analysis and the main text were done by CVB, and contributions were also made by TSu, MH, AL, GWM, UN, CT, IQ, ST, DV, PY, YZ, and TP. TSu, TP, GC, and VA supervised CVB. CT, SSD, and GWM conceptualized the study. The simulations were run by CVB, AL, GWM, UN, GP, IQ, TSe, DV, PY, and YZ. All authors contributed to the discussion of the results and the corrections to the paper.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Atmospheric Chemistry and Physics. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

Disclaimer

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Acknowledgements

We thank Beiping Luo for providing the SAGE-3λ aerosol data used in this work. Furthermore, we thank the editor and the two anonymous referees, who have provided valuable feedback on the article. This work is a contribution to the “Interactive Model Intercomparison Project (ISA–MIP)” from the WCRP/SPARC activity “Stratospheric Sulfur and its Role in Climate (SSiRC)”. ECHAM5–HAM computations and analysis were performed on the computer of the Deutsches Klimarechenzentrum (DKRZ), using resources granted by its Scientific Steering Committee (WLA, under project ID bm0885). NCAR’s Community Earth System Model project has been supported primarily by the National Science Foundation. Computing and data storage resources, including the Cheyenne supercomputer (10.5065/D6RX99HX, National Center for Atmospheric Research, Computational and Information Systems Lab, 2024), were provided by the Computational and Information Systems Laboratory (CISL) at NCAR. SOCOL–AERv2 simulations have been performed at the ETH cluster EULER and Swiss National Supercomputing Centre (CSCS) under project s1191. Valentina Aquila has been supported by the NASA Interdisciplinary Research in Earth Science Program (program manager Ben Phillips) and the NASA Aura Science Team (program manager Richard Eckman). The UM–UKCA simulations were carried out on the UK ARCHER national supercomputing service, with data analysis and storage within the UK collaborative JASMIN data facility. We acknowledge Mohit Dalvi, Nicolas Bellouin, James Brooke, Lauren Marshall, and Kamalika Sengupta for their involvement with the UM–UKCA aerosol code. Slimane Bekki acknowledges the support of CNES (Centre National d’Études Spatiales) to the BeSAFE balloon mission. Earth Simulator was used for the MIROC-CHASER simulations via the support of the Japan Agency for Marine-Earth Science and Technology.

Financial support

Claudia Timmreck and Ulrike Niemeier received support from the Deutsche Forschungsgemeinschaft (DFG) Research Unit VolImpact (FOR2820, grant no. 398006378) within the projects VolARC and VolClim. Yunqian Zhu received funding from NOAA’s Earth Radiation Budget (ERB) Initiative (CPO, grant no. 03-01-07-001) and NOAA (grant nos. NA17OAR4320101 and NA22OAR4320151). Timofei Sukhodolov, Eugene Rozanov, and Christina V. Brodowsky received funding from the Swiss National Science Foundation for this study (grant no. 200020-182239, project POLE (Past, present, and future of the ozone layer)). Timofei Sukhodolov also received support from the Karbacher Fonds, Graubünden, Switzerland. Gabriel Chiodo received funding from the Swiss National Science Foundation via Ambizione (grant no. PZ00P2_180043). Sandro Vattioni has been funded via ETH (grant no. ETH-17 19-2). Graham W. Mann and Sandip S. Dhomse have been funded by NERC highlight topic consortium project SMURPHS (grant no. NE/N006038/1) and the standard grant project MeteorStrat (grant no. NE/R011222/1). Graham W. Mann received funding from the UK National Centre for Atmospheric Science (NCAS) through the NERC multi-center Long-Term Science program on the North Atlantic climate system (ACSIS, NERC grant no. NE/N018001/1) and the Copernicus Atmospheric Monitoring Service (CAMS), one of six services that together form Copernicus, the EU's Earth observation program.

Review statement

This paper was edited by Marc von Hobe and reviewed by two anonymous referees.