Improved tropospheric and stratospheric sulfur cycle in the aerosol–chemistry–climate model SOCOL-AERv2

Abstract

SOCOL-AERv1 was developed as an aerosol–chemistry–climate model to study the stratospheric sulfur cycle and its influence on climate and the ozone layer. It includes a sectional aerosol model that tracks the sulfate particle size distribution in 40 size bins, between 0.39 nm and 3.2 $\mathrm{\mu }\mathrm{m}$. showed that SOCOL-AERv1 successfully matched observable quantities related to stratospheric aerosol. In the meantime, SOCOL-AER has undergone significant improvements and more observational datasets have become available. In producing SOCOL-AERv2 we have implemented several updates to the model: adding interactive deposition schemes, improving the sulfate mass and particle number conservation, and expanding the tropospheric chemistry scheme. We compare the two versions of the model with background stratospheric sulfate aerosol observations, stratospheric aerosol evolution after Pinatubo, and ground-based sulfur deposition networks. SOCOL-AERv2 shows similar levels of agreement as SOCOL-AERv1 with satellite-measured extinctions and in situ optical particle counter (OPC) balloon flights. The volcanically quiescent total stratospheric aerosol burden simulated in SOCOL-AERv2 has increased from 109 Gg of sulfur (S) to 160 Gg S, matching the newly available satellite estimate of 165 Gg S. However, SOCOL-AERv2 simulates too high cross-tropopause transport of tropospheric ${\mathrm{SO}}_{\mathrm{2}}$ and/or sulfate aerosol, leading to an overestimation of lower stratospheric aerosol. Due to the current lack of upper tropospheric ${\mathrm{SO}}_{\mathrm{2}}$ measurements and the neglect of organic aerosol in the model, the lower stratospheric bias of SOCOL-AERv2 was not further improved. Model performance under volcanically perturbed conditions has also undergone some changes, resulting in a slightly shorter volcanic aerosol lifetime after the Pinatubo eruption. With the improved deposition schemes of SOCOL-AERv2, simulated sulfur wet deposition fluxes are within a factor of 2 of measured deposition fluxes at 78 % of the measurement stations globally, an agreement which is on par with previous model intercomparison studies. Because of these improvements, SOCOL-AERv2 will be better suited to studying changes in atmospheric sulfur deposition due to variations in climate and emissions.

1 Introduction

The atmospheric sulfur cycle is of significance for climate, atmospheric chemistry, ecosystems, agriculture, and human health. Sulfate aerosol reflects incoming shortwave solar radiation, leading to a cooling effect at the Earth’s surface. Sulfate aerosol also absorbs outgoing longwave radiation, leading to a warming of the lower stratosphere. In addition to these direct radiative effects, sulfate particles act as cloud condensation nuclei, leading to cloud formation and an indirect radiative effect . Sulfate particles affect air chemistry in the troposphere and stratosphere by catalyzing chemical reactions that deactivate nitrogen . In the cold polar winter stratosphere, they affect ozone depletion by activating chlorine species and serving as condensation nuclei for polar stratospheric clouds . For decades, atmospheric sulfate deposition has been a concern due to its role in acidification of soils and surface waters . On the other hand, sulfur is a macronutrient for plants and decreasing sulfur deposition has led to increased demand for sulfur fertilizers in certain regions . These wide-ranging impacts have motivated the development of atmospheric sulfur models.

Sulfur is emitted to the atmosphere in various compounds through both natural and anthropogenic sources. Natural sources of sulfur include ${\mathrm{SO}}_{\mathrm{2}}$ from eruptive and degassing volcanoes and dimethylsulfide (DMS) from marine phytoplankton and in small amounts from the terrestrial biosphere. Anthropogenic activities such as fossil fuel combustion, metal smelting, and biomass burning release sulfur mainly as ${\mathrm{SO}}_{\mathrm{2}}$ to the atmosphere . Short-lived sulfur compounds such as hydrogen sulfide (${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$), DMS, and carbon disulfide (CS${}_{\mathrm{2}}$) are almost completely oxidized in the troposphere and thus do not enter the stratosphere in large amounts. Longer-lived sulfur compounds, such as carbonyl sulfide (OCS) and to some extent sulfur dioxide (${\mathrm{SO}}_{\mathrm{2}}$), are transported to the stratosphere where they ultimately oxidize to gaseous sulfuric acid (${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$) . ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ molecules nucleate to form new sulfate aerosol particles or condense on existing particles. In the stratosphere, this sustains a layer of binary ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$-${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ particles between 15 and 30 km, often called the “Junge layer” . The Junge layer is intensified compared to background conditions by sporadically occurring volcanic eruptions whose emissions reach the stratosphere. In the troposphere, ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$-containing particles are removed by wet and dry deposition, closing the atmospheric sulfur cycle .

OCS was first suggested to be an important contributor to the stratospheric aerosol layer by . Recent modeling studies have quantified the contribution of different sulfur compounds to the stratospheric aerosol burden. attributed 70 % of the background, nonvolcanic stratospheric aerosol burden above 20 km to OCS oxidation. suggested that 56 % of the stratospheric burden arises due to OCS and 32 % due to ${\mathrm{SO}}_{\mathrm{2}}$ emissions. It must be noted that these studies calculated these contributions by turning off emissions from all other sulfur species, and lower sulfur emissions can lead to smaller aerosol particles, slower sedimentation, and longer aerosol lifetimes.

As reviewed by , other studies have emphasized the minor, yet nonnegligible, contribution of nonsulfate species to the stratospheric aerosol burden. Meteoritic dust particles enter the atmosphere at a rate of 3–300 t d${}^{-\mathrm{1}}$ , which would correspond to 0.15 %–15 % of the stratospheric aerosol mass flux estimated by . Modeling and measurements suggest that organic carbon is a significant fraction of the aerosol mass in the lowest part of the stratosphere. The nonsulfate aerosol species are often not considered in stratospheric modeling studies, despite their possible contribution to observed aerosol quantities.

There are around 15 active models that include stratospheric aerosol microphysics, which can be separated into models with sectional ($\sim \mathrm{1}/\mathrm{3}$) or modal ($\sim \mathrm{2}/\mathrm{3}$) size distributions . SOCOL-AERv1 is a model with a sectional scheme that divides the sulfate aerosol size distribution into 40 bins . The model succeeded in reproducing the observed background stratospheric aerosol extinctions compared to the Stratospheric Aerosol and Gas Experiment II (SAGE II) and Halogen Occultation Experiment (HALOE) measurements , as well as the particle size distributions measured by optical particle counters (OPCs) in the midlatitudes . The SOCOL-AERv1-simulated aerosol burden of 109 Gg sulfur (S) also matched the stratospheric burden calculated from SAGE-$\mathrm{4}\mathit{\lambda }$ data (112 Gg S).

Despite the good agreement of SOCOL-AERv1 with stratospheric aerosol observations, several aspects of the tropospheric sulfur cycle are treated in a coarse manner. For example, the wet and dry deposition schemes are not interactive, i.e., wet removal of precursor gases and aerosol does not depend on the grid cell precipitation and dry deposition does not depend on the land surface type or particle size. The Model Intercomparison Project on the climatic response to volcanic forcing (VolMIP) , which investigated the response of four chemistry–climate models (CCMs) to the 1815 eruption of Mt. Tambora in Indonesia, highlighted several concerns with the deposition fluxes simulated by SOCOL-AERv1. Compared to the other models with interactive deposition schemes, SOCOL-AERv1 displayed lower sulfate deposition in the midlatitude storm tracks since its wet deposition scheme is not linked with precipitation. As well, SOCOL-AERv1 overestimated sulfate deposition to polar ice sheets in both the preindustrial background and Tambora cases . Improvements to the deposition schemes in SOCOL-AER are expected to lead to better reconstructions of past volcanic activity from deposited sulfate.

Since was published, there have also been substantial updates and changes to the stratospheric aerosol observations. OPC measurements have undergone revision due to a correction in the counting efficiency of the instrument . Updated extinction values are available through the Global Space-based Stratospheric Aerosol Climatology (GloSSAC) project . This dataset was used to construct the SAGE-$\mathrm{3}\mathit{\lambda }$ dataset of stratospheric aerosol burdens, which is an input for models in phase 6 of the Coupled Model Intercomparison Project (CMIP6). New in situ measurements of ${\mathrm{SO}}_{\mathrm{2}}$ in the upper troposphere have raised a discussion about the magnitude of the cross-tropopause ${\mathrm{SO}}_{\mathrm{2}}$ flux . Here, we use these newly available datasets to evaluate results from the updated SOCOL-AERv2 model.

This paper outlines the changes that have been made from SOCOL-AERv1 to v2 through the implementation of the new interactive deposition scheme and other improvements. Section  summarizes the changes to the SOCOL-AER code and details the experimental setup of the simulations. Section  discusses the results of three types of simulations: year 2000 time-slice runs (Sect. –), 2000–2010 transient runs (Sect. ), and post-Pinatubo transient runs (Sect. ). We compare the model with stratospheric aerosol observations, from both nonvolcanic background and post-Pinatubo periods, as well as with surface measurements of wet and dry deposition fluxes. Section  presents updated estimates for the year 2000 atmospheric sulfur budget. Section  draws the conclusions of this work.

2 Description of model simulations

2.1 Year 2000 time-slice simulations

The development of SOCOL-AERv2 consisted of corrections to the SOCOL-AERv1 code and implementations of new schemes. To compare directly with the reference simulation from we run time-slice simulations for the year 2000 at each stage of the code changes. For each model run we simulate 10 years, taking the first 5 years as a spin-up from the initial conditions and analyzing only the last 5 years. In this section, we describe changes in the code for each of the time-slice simulations (summarized in Table ).

Table 1

Description of the year 2000 time-slice simulations that were produced during the development of SOCOL-AERv2. Each simulation builds on the previous one, adding new features or correcting a certain process in the model. Two additional sensitivity runs are also listed that reduce the cross-tropopause transport of ${\mathrm{SO}}_{\mathrm{2}}$ and tropospheric sulfate aerosol.

Simulation name	Description	Section
SHENG31	Model version (in T31 resolution) used in	2.1.1
SOCOL-AERv1	SHENG31 model code but using T42 resolution	2.1.1
DRYRAD	Implementation of dry sulfate aerosol radius as the binning scheme instead of wet radius	2.1.2
CONSERVE	Improvement of mass conservation in aerosol microphysical schemes	2.1.3
CCMI	Expansion of tropospheric chemistry scheme	2.1.4
BNDLAYER	Limiting emission and deposition boundary conditions to the lowermost model layer	2.1.5
DRYDEP	Addition of interactive dry deposition scheme	2.1.6
WETDEP	Addition of interactive wet deposition scheme	2.1.7
AQCHEM	Fixes to the aqueous chemistry scheme (time step, transfer to wet deposition flux, cloud pH)	2.1.8
SOCOL-AERv2	Correcting the supercooled liquid fraction to CALIOP-observed values for aqueous chemistry	2.1.9
Additional runs based on SOCOL-AERv2
ICE-OX	Oxidation of S(IV) to S(VI) occurs in cloud ice water in addition to liquid water	2.1.10
AER-SCAV	Aerosol scavenging coefficient on ice is increased by a factor of 20	2.1.10

2.1.1 Rerunning v1 in T31 (simulation: SHENG31) and T42 (SOCOL-AERv1)

A full description of SOCOL-AERv1 can be found in , so here we will only summarize the main aspects of the model. The model originated from the SOCOLv3 chemistry–climate model , which consists of the middle atmosphere version of the global circulation model (GCM) European Centre/Hamburg 5 (MA-ECHAM5) and the chemistry model MEZON . SOCOLv3 includes 39 hybrid vertical levels ranging from the Earth surface up to 0.01 hPa (80 km). The model is an atmosphere-only model prescribing global sea surface temperatures and sea ice coverage with observed data from the Hadley Centre . The quasi-biennial oscillation is produced in the model by relaxing the simulated zonal winds in the equatorial stratosphere to observed wind profiles .

The chemistry module in SOCOL-AERv1 includes a comprehensive range of stratospheric chemical reactions and a simplified set of tropospheric reactions: of the atmospheric hydrocarbons, only methane photochemistry is included. introduced online sulfur chemistry and sulfate aerosol microphysics to the SOCOL model, based on the two-dimensional sulfate aerosol model AER . The model considers eight gaseous sulfur species: OCS, CS${}_{\mathrm{2}}$, ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$, DMS, methanesulfonic acid (MSA), ${\mathrm{SO}}_{\mathrm{2}}$, sulfur trioxide (${\mathrm{SO}}_{\mathrm{3}}$), and ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$. Sulfate aerosol particles are resolved in 40 size bins, ranging in radius from 0.39 nm to 3.2 $\mathrm{\mu }\mathrm{m}$, with sequential bins doubling in volume. Chemical reaction rate coefficients and absorption cross sections of all reactions, including sulfur reactions, follow the recommendations by . In the aqueous phase, sulfur is described as S(IV) and S(VI) without further speciation. Aqueous oxidation of S(IV) by ozone (${\mathrm{O}}_{\mathrm{3}}$) and hydrogen peroxide (${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$) is calculated by the model using the scheme by . The aqueous production flux of S(VI) is added to the atmospheric sulfate aerosol tracers with the flux to each bin proportional to the bin volume.

The microphysical scheme of SOCOL-AERv1 considers the nucleation , composition , growth through ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ condensation, evaporation , coagulation , and sedimentation of sulfate aerosol . SOCOL-AERv1 employs a crude altitude-varying lifetime approach for tropospheric wet removal of sulfur species, with ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ and MSA having a mean column lifetime of 5 d and ${\mathrm{SO}}_{\mathrm{2}}$ having a mean lifetime of 2.5 d . Dry deposition of ${\mathrm{SO}}_{\mathrm{2}}$, MSA, ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$, and sulfate aerosol is modeled assuming a deposition velocity of 1 cm s${}^{-\mathrm{1}}$ at the ground. The model is run with operator splitting so that dynamical quantities are recalculated every 15 min, whereas the chemistry, aerosol microphysics, and radiation schemes are called every 2 h. Twenty sub-time-steps are used for the aerosol microphysical scheme, yielding an aerosol microphysical time step of 6 min.

The model's boundary conditions that we use for the year 2000 time-slice simulations are identical to . ${\mathrm{SO}}_{\mathrm{2}}$ is emitted from anthropogenic and biomass burning sources according to a gridded emission inventory for the year 2000 and from continuous volcanic degassing . DMS fluxes are calculated online using a wind-driven parametrization and a climatology of sea surface DMS concentrations . As in , 1 Tg S yr${}^{-\mathrm{1}}$ of CS${}_{\mathrm{2}}$ is emitted between the latitudes of 52${}^{\circ }$ S and 52${}^{\circ }$ N. The mixing ratios of ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ and OCS are fixed at the surface to 30 pptv and 500 pptv , respectively.

To ensure comparability of results with the new development runs for this paper, we have rerun the source code from in two experiments. As opposed to the simulations in , we output the sulfur cycle burdens and fluxes as accumulated rather than instantaneous quantities to reduce the influence of diurnal cycling on the 12-hourly output of the model. To test the effect of horizontal resolution on the atmospheric sulfur cycle we ran one simulation at T31 resolution ($\sim \mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$ in latitude and longitude, SHENG31) and one simulation at T42 resolution (SOCOL-AERv1).

2.1.2 Dry radius binning scheme (DRYRAD)

In SOCOL-AERv1, the sulfate aerosol is resolved in wet radius bins. Uptake and evaporation of ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ during transport-induced changes in relative humidity and temperature cause shifts in the sulfate mass distribution with respect to wet aerosol radius, the coordinate variable. In SOCOL-AERv1, this process was treated by rescaling the number of sulfate aerosol particles so that each bin would have the correct ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ weight percent. Although this procedure guarantees the conservation of the total number of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ molecules, it does not conserve the aerosol number density. To amend this, SOCOL-AERv2 resolves the sulfate aerosol distribution in dry radius bins, similar to the approach of other sectional models e.g.,. Dry radius bins can also be described as aerosol ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ mass bins.

The new dry radius bins range from 0.39 nm to 3.2 $\mathrm{\mu }\mathrm{m}$, corresponding to 2.8 to $\mathrm{1.6}\times {\mathrm{10}}^{\mathrm{12}}$ molecules of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ per particle (assuming a ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ density of 1.8 g cm${}^{-\mathrm{3}}$), with molecule number doubling between bins. Since aerosol microphysics schemes and heterogeneous chemistry on sulfate aerosol require wet aerosol volume and ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ weight percent, we calculate these quantities for each bin online, based on the grid cell temperature and relative humidity. This change in the dimension variable of SOCOL-AER necessitated several changes to the sulfate condensation and coagulation schemes to ensure that the transfer of aerosol from bin to bin is based on molecular fluxes rather than aerosol volume fluxes. For calculation of aerosol radiative properties, a new look-up table was produced as a function of relative humidity and temperature. To ease interpretation of the output, the outputted aerosol bins of SOCOL-AER are rebinned to the previous wet volume binning approach.

2.1.3 Mass conservation issues (CONSERVE)

In the CONSERVE simulation, corrections were made to the aerosol microphysics scheme in SOCOL-AER, mainly to improve aerosol mass conservation. In the scheme calculating ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ condensation and evaporation, the equation for the “effective” mean free path of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ molecules in air was corrected to agree with Eq. (6) from . An additional constraint was added in the ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ condensation scheme such that the flux of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ from the gas phase must equal the flux of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ into the particle phase. This improves mass conservation in cases when ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ is depleted below the saturation vapor pressure within one time step. Furthermore, the aerosol sedimentation scheme in SOCOL-AERv1 was not applied within boundary layer levels and sedimenting aerosol from the model level above the boundary layer was artificially removed. In the CONSERVE simulation this is amended: sedimenting aerosol from the model level above the boundary layer is added to the layer below. had implemented several forced mass conservation checks on the total (gas-phase and aerosol) ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ burden in each grid cell. If the total burden had changed by more than 0.1 % during aerosol microphysics, ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ aerosol and gas-phase mixing ratios would be scaled to agree with the total ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ burden before microphysics calculations. These forced mass corrections were found to be unnecessary after the above improvements to the microphysics scheme, and therefore they were removed from SOCOL-AERv2.

2.1.4 Merging CCMI additions with SOCOL-AER (CCMI)

Since the publication of , improvements have been made to the SOCOL model in preparation of the coordinated simulations within the Chemistry Climate Model Initiative (CCMI), mainly related to the improvement in tropospheric chemistry processes . Many of these improvements have been merged into SOCOL-AERv2 and have upgraded the representation of chemistry in our model, in particular in the troposphere.

In SOCOL-AERv1, as well as SOCOLv3 , ozone-depleting substances (ODSs) were transported in three families (short-lived Cl, long-lived Cl, and Br) to save computational cost. With modern supercomputers this treatment is no longer necessary and the ODS species are transported individually. The individual treatment of ODS species avoids a repartitioning of the family members, based on simplified age-of-air estimates, after each transport step. The chemistry scheme was expanded in the CCMI simulation, as described in . We included the Mainz Isoprene Mechanism (MIM-1) in SOCOL-AER. This scheme considers the degradation of isoprene and necessitates the addition of 14 organic species and 44 chemical reactions to SOCOL-AER . Additional CO emissions were added to the model to account for the effect of oxidation of nonmethane volatile organic compounds (NMVOC) emitted from anthropogenic, biogenic, and biomass burning sources. Lightning NO${}_x$ is now calculated interactively based on cloud-top height and grid cell scaling factors from satellite observations . A cloud modification factor approach was implemented to account for the effect of clouds on photolysis rates. We derived a new look-up table of photolysis rates averaged over two solar cycles (22 years) from a comprehensive reconstruction of total and spectral solar irradiance (NRLSSI) by , which was used in the CCMI REF-C1 experiment. Additional heating through Hartley and Huggins bands of ozone has also been implemented into SOCOL-AER. As documented by , ${\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$ hydrolysis on tropospheric aerosols is now included in SOCOL-AER. Methanesulfonic acid (MSA) chemistry is solved in the explicit scheme instead of the implicit Newton–Raphson scheme since otherwise the chemical solver does not properly converge. Reaction rates have been updated or added from the NASA JPL data evaluation no. 18 .

2.1.5 Treatment of the boundary layer (BNDLAYER)

In the previous simulations, to allow for rapid boundary layer mixing, emissions of chemical species were immediately dispersed over the four lowermost model levels ($\sim \mathrm{1}$ km altitude), species with prescribed mixing ratios (including OCS and ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$) were dispersed over the six lowermost levels ($\sim \mathrm{2}$ km altitude), and dry deposition of species occurred out of the four lowermost model levels. While such a coarse approach was sufficient for stratospheric applications, it is inadequate for deposition flux and tropospheric lifetime calculations. Instead of emitting the species in multiple lower layers, SOCOL-AERv2 emits only in the first model layer ($\sim \mathrm{70}$ m), from which the species are mixed via the model's boundary layer parametrizations. The BNDLAYER simulation tests specifically the effect of using only one model layer for emission and prepares the model for revised dry deposition boundary conditions.

2.1.6 Interactive dry deposition (DRYDEP)

We implemented the interactive dry deposition scheme described in in SOCOL-AERv2, replacing the simple prescribed constant deposition velocities of SOCOL-AERv1. The new treatment is based on the DRYDEP scheme in the EMAC model . Dry deposition velocities are calculated using an interactive resistance-based approach which considers surface properties, the solubility and reactivity of each gas tracer, and the radius and density of aerosol tracers . Effective Henry’s law constants for near-neutral pH and reactivity of gas tracers are taken from . These improvements are tested in the DRYDEP simulation.

2.1.7 Interactive wet deposition (WETDEP)

An interactive wet deposition scheme was added to SOCOL-AER, based on the SCAV submodule in the EMAC model . Grid scale variables from ECHAM5 such as liquid and ice water contents, cloud cover, convective and large-scale rain and ice formation and precipitation fluxes, and the convective upward mass flux are used by the wet deposition scheme. Since our model does not include a comprehensive cloud aqueous chemistry mechanism, we implemented the EASY2 version of the SCAV submodule with a constant pH of 5 for cloud water and rain water. The constant pH of 5 is within the wide range of pH values (3.6–7) measured by several hill cap cloud field campaigns . Scavenging coefficients for gas-phase species are calculated based on Henry’s law equilibrium constants. Scavenging of aerosol is based on a radius-dependent calculation of nucleation and impaction scavenging. During cloud evaporation, all scavenged gas-phase species are released to the atmosphere in their original species, whereas evaporating scavenged sulfate aerosol species are transferred to the largest aerosol size bin. The wet deposition scheme is applied to ${\mathrm{SO}}_{\mathrm{2}}$, gaseous ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$, and sulfate aerosol, as well as other gas chemical tracers such as ${\mathrm{O}}_{\mathrm{3}}$, ${\mathrm{HNO}}_{\mathrm{3}}$, ${\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$, ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$, etc. The WETDEP simulation includes this new interactive wet deposition scheme instead of the fixed wet deposition lifetimes used in SOCOL-AERv1.

2.1.8 Improvement in aqueous-phase chemistry (AQCHEM)

In the ${\mathrm{SO}}_{\mathrm{2}}$ aqueous chemistry subroutine of SOCOL-AERv1, the pH of clouds is prescribed vertically according to so that pH equals 3 from the surface to 600 hPa and 4.5 above 600 hPa. However, this paper reported modeled pH within a single cumulus cloud, where liquid water content increased with height. This pH distribution is therefore not applicable to the whole atmosphere. Although in the interactive wet deposition scheme a constant cloud pH of 5 is used (Sect. ), aqueous-phase sulfur chemistry is more sensitive to the choice of pH and therefore a more detailed pH distribution was applied. We use an approximation of the modeled cloud pH from for the revised aqueous chemistry routine. Between the surface and 600 hPa, north of 20${}^{\circ }$ N a cloud pH of 5.2 is used and south of 20${}^{\circ }$ N a cloud pH of 4.2 is used. Above the 600 hPa level, a uniform pH of 3.5 is used.

In SOCOL-AERv1, the ${\mathrm{SO}}_{\mathrm{2}}$ oxidized in the aqueous cloud phase is released as aerosol. With the new interactive wet deposition scheme, it is possible to transfer the oxidized ${\mathrm{SO}}_{\mathrm{2}}$ directly to the scavenged aerosol flux in cloud water. The wet deposition routine is called at each dynamical time step (15 min), while the aqueous-phase chemistry was called at each chemical time step in SOCOL-AERv1 (2 h). To transfer the oxidized ${\mathrm{SO}}_{\mathrm{2}}$ flux to the wet deposition scheme directly, it was both logical and technically simpler to synchronize these two processes. In the AQCHEM simulation, the above changes were added to SOCOL-AER and the aqueous-phase chemistry is called at each dynamical time step.

2.1.9 Final development run for SOCOL-AERv2

Section – complete the description of improvements in developing SOCOL-AERv2 with one exception, namely how ${\mathrm{SO}}_{\mathrm{2}}$ oxidation is calculated in clouds in the mixed-phase temperature regime. This is important because the S(IV) to S(VI) conversion only occurs in the liquid phase in the model. Therefore, in the final development simulation (SOCOL-AERv2) we wanted to investigate whether the aqueous ${\mathrm{SO}}_{\mathrm{2}}$ reaction was hampered by the model's representation of the liquid fraction in mixed-phase clouds. It has recently been discussed in the literature that many general circulation models (GCMs) underpredict the supercooled liquid fraction (SLF) observed by satellite products . The modeled liquid fraction in SOCOL-AERv1 underestimates the fitted CALIOP satellite measurements from throughout most of the mixed-phase cloud temperature range (Fig. S1 in the Supplement). In the SOCOL-AERv2 run, we correct for the influence of the underestimated supercooled liquid fraction on ${\mathrm{SO}}_{\mathrm{2}}$ aqueous chemistry. The ECHAM5 calculated liquid water content (LWC) and ice water content (IWC) are added together and inputted into the aqueous-phase chemistry subroutine as total water content (TWC). If the grid cell temperature ($T$) is in the mixed-phase cloud regime ($-\mathrm{38}$ to 0 ${}^{\circ }$C), we calculate the observed supercooled liquid fraction (SLF) from the fitted sigmoid function from : $$\begin{array}{}\text{1}& {\displaystyle }{\text{SLF}}_{\mathrm{Hu}}={\left[\mathrm{1}+\exp \left(-f\left(T\right)\right)\right]}^{-\mathrm{1}},\text{2}& {\displaystyle }\begin{array}{rl}& f\left(T\right)=\mathrm{5.3608}+\mathrm{0.4025}T+\mathrm{0.08387}{T}^{\mathrm{2}}+\mathrm{0.007182}{T}^{\mathrm{3}}\\ & +\mathrm{2.39}\times {\mathrm{10}}^{-\mathrm{4}}{T}^{\mathrm{4}}+\mathrm{2.87}\times {\mathrm{10}}^{-\mathrm{6}}{T}^{\mathrm{5}},\end{array}\end{array}$$ where $T$ is the air temperature in degrees Celsius (${}^{\circ }$C). The determined SLF can is used to correct LWC in the aqueous chemistry subroutine, i.e., $\text{LWC}={\text{SLF}}_{\text{Hu}}\times \text{TWC}$.

2.1.10 Additional sensitivity runs (ICE-OX, AER-SCAV)

We ran two additional simulations to probe whether the remaining disagreement between observations and SOCOL-AERv2 could be caused by overestimated cross-tropopause fluxes of ${\mathrm{SO}}_{\mathrm{2}}$ and sulfate aerosol. In ICE-OX, the aqueous-phase oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ was allowed to occur in ice water as well as liquid water. Increased oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ in the upper troposphere reduces its cross-tropopause flux. In AER-SCAV, the scavenging coefficient of aerosol particles on ice clouds was increased by a factor of 20 from 0.05 to 1. This enhances the removal of sulfate aerosol in the upper troposphere.

2.2 Years 2000–2010 transient simulations

In order to compare simulated deposition with observations, the model codes from SOCOL-AERv1 and v2 were used to run two sets of transient simulations from 2000 to 2010. Five ensemble members were simulated for both versions of SOCOL-AER and plotted results show ensemble means and standard deviations. For the transient simulations we made several updates to the boundary conditions used in . Anthropogenic emissions were taken from the Community Emissions Data Systems (CEDS), which will be used for CMIP6 simulations . updated the marine DMS dataset to include three times as many DMS measurements as the previous dataset used by . Transient degassing volcanic ${\mathrm{SO}}_{\mathrm{2}}$ emissions were taken from . To represent eruptive emissions, we applied a satellite-derived dataset from . The other data sources for the boundary conditions remained the same as in the time-slice simulations; however, transient boundary conditions were included rather than applying repeating year 2000 values.

2.3 Pinatubo transient simulations

To verify the updated model's performance under volcanically perturbed conditions, we have repeated two experiments from , modeling the Mt. Pinatubo eruption with 7 and 6 Tg S emitted as ${\mathrm{SO}}_{\mathrm{2}}$. As in , we simulated five ensemble members with sulfur mass released from 14 to 15 June 1991 and spread between 16 and 30 km. We performed two additional runs with SOCOL-AERv1 and SOCOL-AERv2, including the Pinatubo eruption magnitude of 7 Tg S but with all other sulfur sources switched off, to check the sulfur mass conservation by analyzing the integrated deposition fluxes. To compare with modeled burdens, we used observational estimates from SAGE-$\mathrm{4}\mathit{\lambda }$ and SAGE-$\mathrm{3}\mathit{\lambda }$ datasets and from the High-Resolution Infrared Radiation Sounder (HIRS) measurements .

3 Results and discussion

3.1 Impacts of performed changes in the development of SOCOL-AERv2

In the following section, we discuss the relevant impacts of each stage of code changes on the atmospheric sulfur cycle. Table  lists the stratospheric and tropospheric burdens of ${\mathrm{SO}}_{\mathrm{2}}$ and sulfate aerosol and total deposition fluxes for each time-slice simulation.

Table 2

Global annual mean sulfate aerosol and ${\mathrm{SO}}_{\mathrm{2}}$ burdens in the troposphere and stratosphere for all volcanically quiescent time-slice (year 2000) simulations in the development of SOCOL-AERv2. Total sulfur deposition (last column) is listed as a check of whether the mass balance of the model has changed. Model results are compared to two datasets, SAGE-$\mathrm{3}\mathit{\lambda }$ and SAGE-$\mathrm{4}\mathit{\lambda }$, which derived stratospheric aerosol burdens from satellite extinction measurements. The model tropopause is used to separate tropospheric and stratospheric burdens, whereas for the SAGE-$\mathrm{3}\mathit{\lambda }$ and SAGE-$\mathrm{4}\mathit{\lambda }$ calculations the MERRA tropopause is used. The WMO assessment calculates a multimodel mean total sulfur deposition flux of $\mathrm{84.8}\pm \mathrm{6.1}$ Tg yr${}^{-\mathrm{1}}$ ($\pm \mathit{\sigma }$) and total emissions of $\mathrm{91.0}\pm \mathrm{7.3}$ Tg yr${}^{-\mathrm{1}}$ for 2001. Simulated interannual variability in one simulation is on the order of 1 %–2 % for burdens and 0.3 % for total deposition.

Simulation	Tropospheric	Stratospheric	Tropospheric	Stratospheric	Total sulfur
	aerosol burden	aerosol burden	${\mathrm{SO}}_{\mathrm{2}}$ burden	${\mathrm{SO}}_{\mathrm{2}}$ burden	deposition
	(Gg S)	(Gg S)	(Gg S)	(Gg S)	(Tg S yr${}^{-\mathrm{1}}$)
SHENG31	395	114	259	12.4	99.8
SOCOL-AERv1	397	116	261	11.9	99.8
DRYRAD	367	125	261	12.0	97.9
CONSERVE	382	128	261	12.3	101.1
CCMI	392	130	242	7.1	101.0
BNDLAYER	236	106	137	5.8	92.5
DRYDEP	508	128	218	6.6	92.8
WETDEP	769	202	351	8.2	93.6
AQCHEM	667	166	245	6.6	94.5
SOCOL-AERv2	640	160	217	6.3	94.4
Additional sensitivity tests
ICE-OX	613	92	188	3.9	94.2
AER-SCAV	579	133	215	6.4	94.5
Observational datasets
SAGE-$\mathrm{4}\mathit{\lambda }$		$\mathrm{117}\pm \mathrm{8}$
SAGE-3$\mathit{\lambda }$		$\mathrm{165}\pm \mathrm{11}$
WMO assessment					$\mathrm{84.8}\pm \mathrm{6.1}$

3.1.1 Rerunning SOCOL-AERv1 in T31 and T42 resolutions

Since the SHENG31 simulation was rerun for this study, several quantities differ slightly in SHENG31 compared to (e.g., 114 vs. 109 Gg S of stratospheric sulfate aerosol). The differences in the quantities could be caused by the switch in output format from instantaneous 12-hourly values to mean 12-hourly values. However, the changes are minor and the overall picture for the sulfur cycle remains unchanged. Refining the horizontal resolution to T42 in SOCOL-AERv1 does not result in substantial changes for the tropospheric and stratospheric aerosol burdens.

3.1.2 Dry radius binning (DRYRAD)

The change from wet radius to dry radius binning reduces the tropospheric aerosol burden and increases the stratospheric aerosol burden in DRYRAD (Table ). The increased stratospheric aerosol burden can be explained by a decrease in the effective aerosol radius, leading to a longer stratospheric lifetime. The decrease in tropospheric aerosol burden occurs mainly around the Northern Hemisphere midlatitudes, where anthropogenic emissions of ${\mathrm{SO}}_{\mathrm{2}}$ are high and thus sulfate particles can grow through condensation. In the DRYRAD version of the model the wet particle radius is no longer restricted to 3.2 $\mathrm{\mu }\mathrm{m}$; accounting for the uptake of water the maximum radius can reach well above 10 $\mathrm{\mu }\mathrm{m}$. The possibility of larger particle formation can lead to enhanced sedimentation velocities and therefore reduced aerosol lifetimes.

3.1.3 Mass conservation fixes (CONSERVE)

The corrections in the sedimentation and ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ condensation schemes do improve the mass conservation of sulfur species. The tropospheric sulfate aerosol burden increases by 4 %, mainly due to the correction of the artificial removal of particles sedimenting from the model level above the boundary layer level. This artificial loss due to sedimentation represents a $\sim \mathrm{3}$ Tg S yr${}^{-\mathrm{1}}$ sink since the outputted total sulfur deposition increases by this amount from the DRYRAD to the CONSERVE simulations. Furthermore, in the DRYRAD simulation the total sum of tropospheric aerosol influxes and outfluxes result in an imbalance of 3339 Gg S yr${}^{-\mathrm{1}}$, which corresponds to about 8 % of the source flux of tropospheric sulfate aerosol. In the CONSERVE simulation, this imbalance is reduced to 63 Gg S yr${}^{-\mathrm{1}}$, i.e., around 0.1 % of the aerosol source flux. These improvements to the model provide more confidence to the outputted sulfur cycle fluxes, which will be used to study the sulfur budget in Sect. .

3.1.4 CCMI chemistry changes (CCMI)

Including the expanded chemistry set and updates from the CCMI version of SOCOLv3 in SOCOL-AER leads to altered distributions of gas-phase species in the troposphere. The relevant change for the tropospheric sulfur cycle is increased mixing ratios of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$, causing increased aqueous conversion of ${\mathrm{SO}}_{\mathrm{2}}$ to S(VI). For this reason, the CCMI simulation shows a 19 Gg S lower ${\mathrm{SO}}_{\mathrm{2}}$ burden and a 10 Gg S larger sulfate aerosol burden in the troposphere than the CONSERVE simulation. Larger OH mixing ratios in the upper troposphere and lower stratosphere (UTLS) reduce the ${\mathrm{SO}}_{\mathrm{2}}$ lifetime, causing 5 Gg lower ${\mathrm{SO}}_{\mathrm{2}}$ and 2 Gg higher sulfate aerosol burdens in the stratosphere. These chemical changes also lead to differences in the Pinatubo simulation, to be discussed in Sect. .

3.1.5 Boundary layer levels (BNDLAYER)

In BNDLAYER the boundary layer conditions are only implemented for the lowest level of the model. The confinement of the boundary layer conditions to one model level reduces the burdens of ${\mathrm{SO}}_{\mathrm{2}}$ and sulfate aerosol in the troposphere and stratosphere. This is a strong effect with reductions of the tropospheric aerosol burden by 40 % and the stratospheric aerosol burden by 20 % in BNDLAYER. The first cause is the reduction in effective S emissions into the atmosphere since ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$ is prescribed to be 30 ppt in only one model level instead of six model levels. Assuming steady-state conditions over the 5-year averaging period, the 8.5 Tg S yr${}^{-\mathrm{1}}$ decrease in total sulfur deposition (Table ) corresponds to an 8.5 Tg S yr${}^{-\mathrm{1}}$ decrease in the sulfur emissions. Another cause for the ${\mathrm{SO}}_{\mathrm{2}}$ and aerosol burden decrease is that ${\mathrm{SO}}_{\mathrm{2}}$ is emitted close to the surface in BNDLAYER, leading to less dispersion of ${\mathrm{SO}}_{\mathrm{2}}$ in the atmosphere and enhanced dry deposition close to emission regions. The shorter ${\mathrm{SO}}_{\mathrm{2}}$ lifetime reduces its atmospheric burden, as well as reducing the conversion of ${\mathrm{SO}}_{\mathrm{2}}$ to ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ and subsequent sulfate aerosol formation.

The correct treatment of the lowermost model levels remains difficult and is model dependent. Owing to their coarse resolution, CCMs cannot resolve the transport of chemical species by rapid boundary layer convection and turbulence. This leaves the boundary layer parametrizations in SOCOL-AER imperfect and the number of model levels that should be included in the emission boundary conditions uncertain. For the subsequent simulations we use the single-level boundary layer treatment.

Figure 1

Annual mean dry deposition velocities for ${\mathrm{SO}}_{\mathrm{2}}$ (a) and sulfate aerosol (b) simulated by SOCOL-AERv2. The mean dry deposition velocity for aerosol particles is calculated by weighting the dry deposition velocity for each size bin with the bin’s mass concentration at the surface. The color bar highlights differences between the newly simulated deposition velocities and the former homogeneous deposition velocity of 1 cm s${}^{-\mathrm{1}}$, which is shown in whitish colors on both plots.

[Figure omitted. See PDF]

3.1.6 Interactive dry deposition (DRYDEP)

The implementation of interactively calculated dry deposition velocities, compared to the previously included constant dry deposition velocities, results in much longer dry deposition lifetimes for both ${\mathrm{SO}}_{\mathrm{2}}$ and sulfate aerosol. ${\mathrm{SO}}_{\mathrm{2}}$ dry deposition velocities decrease more drastically over land than over ocean in the DRYDEP simulation (Fig. ). Over land, the ${\mathrm{SO}}_{\mathrm{2}}$ dry deposition velocity is smaller than 1 cm s${}^{-\mathrm{1}}$, which was the original value set in SOCOL-AERv1. The only locations where ${\mathrm{SO}}_{\mathrm{2}}$ dry deposition velocities are greater than 1 cm s${}^{-\mathrm{1}}$ are above certain parts of the ocean, due to the high solubility of ${\mathrm{SO}}_{\mathrm{2}}$ in waters at near-neutral pH. The dry deposition velocities of ${\mathrm{SO}}_{\mathrm{2}}$ agree well with the distribution simulated by the EMAC model . The resultant longer ${\mathrm{SO}}_{\mathrm{2}}$ dry deposition lifetime increases the tropospheric ${\mathrm{SO}}_{\mathrm{2}}$ burden, the conversion of ${\mathrm{SO}}_{\mathrm{2}}$ to aerosol, and consequently the tropospheric aerosol burden. In addition, dry deposition velocities of sulfate aerosol decrease globally compared to the assumed constant deposition velocity in SOCOL-AERv1 (1 cm s${}^{-\mathrm{1}}$), leading to longer aerosol lifetimes with respect to dry deposition. Due to augmented transport of tropospheric ${\mathrm{SO}}_{\mathrm{2}}$ and primary sulfate aerosol from the troposphere, the stratospheric aerosol burden increases by 22 to 128 Gg S. The changes in DRYDEP largely compensate the changes in BNDLAYER, for which emissions were confined to a single model level.

3.1.7 Interactive wet deposition (WETDEP)

When the constant wet deposition lifetimes for sulfur species are replaced with interactively calculated wet removal in the WETDEP simulation, the ${\mathrm{SO}}_{\mathrm{2}}$ wet deposition flux is reduced from 20.1 to 0.3 Tg yr${}^{-\mathrm{1}}$, revealing an overestimation in the approach of SOCOL-AERv1. With the elimination of the wet deposition sink for ${\mathrm{SO}}_{\mathrm{2}}$, the tropospheric ${\mathrm{SO}}_{\mathrm{2}}$ burden increases by around 60 % and the total (aqueous $+$ gas phase) conversion flux of ${\mathrm{SO}}_{\mathrm{2}}$ to S(VI) increases by around 40 %. In , the mean wet deposition lifetime for ${\mathrm{SO}}_{\mathrm{2}}$ was selected as 2.5 d following the two-dimensional AER model. However, the AER model includes the 2.5 d lifetime for ${\mathrm{SO}}_{\mathrm{2}}$ to account for aqueous oxidation of ${\mathrm{SO}}_{\mathrm{2}}$, which is not explicitly modeled by AER . As SOCOL-AERv1 already includes a mechanism for aqueous oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ by ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$ and ${\mathrm{O}}_{\mathrm{3}}$ in clouds, this resulted in double counting of the loss of ${\mathrm{SO}}_{\mathrm{2}}$ by aqueous oxidation in previous simulations. The aerosol wet deposition maps and the relative difference between the DRYDEP and WETDEP simulations are shown in Fig. . The inclusion of an interactive wet deposition enhances sulfate aerosol deposition in areas with high precipitation and suppresses it in drier regions. Sulfate deposition is reduced over polar regions, the eastern part of ocean basins, and the Sahara (lower precipitation regions), and is enhanced in the tropics and midlatitude storm tracks (higher precipitation regions). The reductions in sulfate deposition fluxes above Greenland and Antarctica are notable since SOCOL-AERv1 overestimated the magnitude of polar sulfate deposition fluxes compared to ice core measurements, which are used as proxies for past volcanic eruptions . Calculating a global aerosol wet deposition lifetime with respect to wet deposition (lifetime is tropospheric aerosol burden divided by aerosol wet deposition flux), DRYDEP has an aerosol wet deposition lifetime of 4.9 d and WETDEP has a lifetime of 5.1 d. Therefore, there is not a large change in the global aerosol wet deposition lifetime; however, the spatial distribution of the wet deposition sink has shifted. In WETDEP the tropospheric column wet deposition lifetime of sulfate aerosol varies from 2 d in the northern midlatitude storm tracks to more than 3 years over the southwestern United States. The introduction of interactive wet deposition to SOCOL-AER has the largest impact of any step on the stratospheric sulfate burden. Driven by the longer ${\mathrm{SO}}_{\mathrm{2}}$ wet deposition lifetime, the stratospheric sulfate aerosol burden climbs by around 60 % to 202 Gg S. As will be discussed in Sect. , this value is much higher than the inferred stratospheric burden from SAGE II data for background nonvolcanic conditions. To improve the agreement with observations, we focus on a possible underestimation of the sulfate aqueous chemistry flux since the unintended double counting of this flux led to good agreement of SOCOL-AERv1 with stratospheric observations .

Figure 2

Maps of sulfate wet deposition simulated by SOCOL-AER for the year 2000. Three plots are shown: wet deposition in the DRYDEP simulation, which uses the old wet deposition scheme (a); wet deposition in the WETDEP simulation, which uses the interactive wet deposition scheme (b); and relative percent differences in deposition, which is WETDEP minus DRYDEP (c).

[Figure omitted. See PDF]

3.1.8 Aqueous chemistry changes (AQCHEM)

In the AQCHEM simulation we amended the cloud pH distribution for aqueous chemistry, reduced the aqueous chemistry time step to 15 min, and directly transferred oxidized S(IV) to the scavenged sulfate aerosol in the wet deposition scheme. This increases the aqueous oxidation flux of S(IV) to S(VI) by around 50 %. The enhanced aqueous conversion of ${\mathrm{SO}}_{\mathrm{2}}$ to sulfate aerosol leads to increased aerosol formation in the lowermost troposphere, where deposition is efficient. This results in smaller tropospheric burdens of both ${\mathrm{SO}}_{\mathrm{2}}$ and sulfate aerosol, meaning that there is also less transport of tropospheric S to the stratosphere. The stratospheric aerosol burden decreases by 18 % from 202 to 165 Gg S. From separate sensitivity studies (not shown), we find that the shorter aqueous chemistry time step is the main cause of the increased aqueous flux. investigated the sensitivity of sulfate aqueous chemistry to different settings and also found that reducing the aqueous chemistry time step increases the conversion of S(IV) to S(VI). This is because the Henry’s law equilibration rate and aqueous oxidation is fast so with shorter time steps more ${\mathrm{SO}}_{\mathrm{2}}$ can be dissolved in cloud droplets and converted to S(VI).

3.1.9 SOCOL-AERv2 and aqueous chemistry in the supercooled liquid fraction

Because of the underprediction of the SLF (Fig. S1) and oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ occurring only in liquid water in SOCOL-AER, the oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ is likely underestimated in the upper troposphere, leading to a too intensive transport of ${\mathrm{SO}}_{\mathrm{2}}$ to the stratosphere. Therefore, in the SOCOL-AERv2 simulation we increased the supercooled liquid fraction to agree with the SLF–temperature relationship observed by CALIOP . This increase in SLF enhances the ${\mathrm{SO}}_{\mathrm{2}}$ oxidation rate in the middle and upper tropospheric mixed-phase clouds, reducing the cross-tropopause ${\mathrm{SO}}_{\mathrm{2}}$ flux by around 10 %. However, the impact on the stratospheric aerosol burden is minor, with only a reduction of 6 Gg S ($-\mathrm{4}$ %) compared to the AQCHEM simulation. Therefore, the underestimation of the SLF does not play a major role in SOCOL-AER's stratospheric sulfur cycle. However, the amount of ${\mathrm{SO}}_{\mathrm{2}}$ oxidation in the upper troposphere may be affected by other processes, e.g. by oxidation on ice surfaces. This will be discussed further in Sect. .

3.2 Comparison of SOCOL-AER versions with stratospheric observations

In this section, we will compare the SOCOL-AERv1 and v2 simulations with observations to understand how the model results change in the new version and where deficiencies remain.

3.2.1 Comparison with SAGE II-derived burdens

compared the modeled stratospheric sulfate aerosol burden to the value calculated by the SAGE-$\mathrm{4}\mathit{\lambda }$ method . In this method, extinctions measured by the SAGE II satellite product are used to estimate the stratospheric aerosol size distribution, which can then be used to determine the aerosol burden. The background stratospheric aerosol burden derived from SAGE-$\mathrm{4}\mathit{\lambda }$ almost exactly matched the burden simulated by SOCOL-AERv1. Since that time, a new SAGE II retrieval has been published as part of the GloSSAC database . A new method (SAGE-$\mathrm{3}\mathit{\lambda }$) has been used to calculate the aerosol size distribution from the GloSSAC database. In this method, the surface area density and mass density of very small particles, which are invisible to the satellite extinction measurements, are added to the lognormal size distributions derived from the GloSSAC data. The stratospheric aerosol burden derived from SAGE-$\mathrm{3}\mathit{\lambda }$ for the volcanically quiescent period 2000–2004 is $\mathrm{165}\pm \mathrm{11}$ ($\pm \mathit{\sigma }$) Gg S. This aerosol burden is about 40 % larger than the stratospheric burden calculated from SAGE-$\mathrm{4}\mathit{\lambda }$, $\mathrm{117}\pm \mathrm{8}$ Gg S. (Note that this value differs slightly from the value reported in , 112.5 Gg S, possibly due to different assumptions about the tropopause height.) The addition of small aerosol particles derived from OPC measurements contributes 18 Gg S to the SAGE-$\mathrm{3}\mathit{\lambda }$. The rest of the increase from SAGE-$\mathrm{4}\mathit{\lambda }$ to SAGE-$\mathrm{3}\mathit{\lambda }$ can be attributed to the new retrieval methods.

The stratospheric aerosol burden simulated by SOCOL-AERv2, 160 Gg S, agrees well with the SAGE-$\mathrm{3}\mathit{\lambda }$-derived burden of 165 Gg S. However, evaluating a model's performance with the stratospheric aerosol burden is not straightforward since both SAGE-$\mathrm{3}\mathit{\lambda }$ and SAGE-$\mathrm{4}\mathit{\lambda }$ are themselves derived products and not direct measurements. The retrieval of size distributions from measured SAGE II wavelengths is uncertain, as can be seen when comparing the change between SAGE-$\mathrm{3}\mathit{\lambda }$ and SAGE-$\mathrm{4}\mathit{\lambda }$ stratospheric burdens. In addition, the MERRA climatological tropopause height was used to calculate the stratospheric burden for the SAGE-$\mathrm{3}\mathit{\lambda }$ and SAGE-$\mathrm{4}\mathit{\lambda }$ products. However, for SOCOL-AER's burden, the WMO-defined (World Meteorological Organization) tropopause height from the model was used . Differences between the tropopause heights used in different calculations can play a big role in the derived burden since the majority of the aerosol burden is located in the lower stratosphere. For example, if the tropopause height from SOCOL-AER instead of MERRA is used, the stratospheric burdens derived from SAGE-$\mathrm{3}\mathit{\lambda }$ and SAGE-$\mathrm{4}\mathit{\lambda }$ are around 7 % smaller. For these reasons we will evaluate SOCOL-AER with the extinctions measured directly by SAGE II of version GloSSACv1.0, in addition to the derived burdens.

3.2.2 Comparison with SAGE II extinctions

Figure  shows the comparison of annual mean SOCOL-AERv1 and v2 extinctions with SAGE II at the Equator and 45${}^{\circ }$ N for 525 and 1020 nm. Below 20 km at the Equator, SOCOL-AERv2 shows higher extinctions at both wavelengths, which match better with the SAGE II observations. However, in the lowest 1–3 km of the stratosphere, organics are a nonnegligible fraction of the overall aerosol burden and therefore can contribute to the aerosol extinction observed by SAGE II. If anything, SOCOL-AER as a sulfate aerosol-only model should underestimate the extinction at these altitudes, although to what extent is unknown. Since SOCOL-AERv2 matches or overestimates the SAGE II extinctions in the lowermost stratosphere, SOCOL-AERv2 may have too high sulfate aerosol concentrations in the lower stratosphere. Between 20 and 25 km SOCOL-AERv2 overestimates the SAGE II extinctions, while SOCOL-AERv1 matches observations. Between 25 and 30 km both model versions overestimate the SAGE II extinctions. The model versions are within observed variability between 30 and 35 km; however, above 35 km they tend to underestimate the extinction, possibly because of meteoritic dust, which is the major contributor to extinction in the upper stratosphere . The comparison at 45${}^{\circ }$ N is similar to the equatorial comparison with SOCOL-AERv2 overestimating the observed aerosol extinctions below 20 km and otherwise showing similar behavior to SOCOL-AERv1.

Since the overestimation of SOCOL-AERv2 in the lower stratosphere originates from the introduction of interactive deposition schemes, it possibly stems from too fast cross-tropopause transport of primary sulfate aerosol and/or ${\mathrm{SO}}_{\mathrm{2}}$, whereas the better agreement of SOCOL-AERv1 may be fortuitous due to the double counting of the ${\mathrm{SO}}_{\mathrm{2}}$ oxidation flux in the wet-deposition scheme. SOCOL-AERv2 is the version that is more physically consistent in its representation of the tropospheric sulfur cycle. However, several outstanding issues remain in SOCOL-AERv2's representation of sulfate aerosol extinction below 20 km at 45${}^{\circ }$ N and between 20 and 30 km at the Equator.

Figure 3

Comparison between annual mean model extinctions at 525 and 1020 nm and SAGE II measurements from the GloSSAC project at the Equator (a) and 45${}^{\circ }$ N (b). Observations are averaged between 2000 and 2004, representing the volcanically quiescent part of the record. Model results are averaged over 5 years of the year 2000 time slice for SOCOL-AERv1 and SOCOL-AERv2. Horizontal bars represent the modeled or observed standard deviation. The highlighted region in the upper plot corresponds to the altitudes where nonsulfate aerosols may play a role.

[Figure omitted. See PDF]

3.2.3 Comparison with OPC size distributions

We also compare the SOCOL-AER simulations with in situ OPC measurements from Laramie, USA, and Lauder, New Zealand (Fig. ). Since the publication of , the counting efficiencies of OPC channels as a function of radius have undergone important revisions . In Fig. , we apply the measured counting efficiencies for the channels $r>\mathrm{0.15}$, $r>\mathrm{0.25}$, and $r>\mathrm{0.30}$ $\mathrm{\mu }\mathrm{m}$ from to the SOCOL-AER size bins (counting efficiencies were not measured for other channels). In this manner, we can calculate the number density that an OPC instrument would measure given a simulated size distribution.

Figure 4

Number densities of particle size bins measured by OPC and modeled by SOCOL-AERv1 and v2 over Laramie, Wyoming, USA (41${}^{\circ }$ N, 105${}^{\circ }$ W), and Lauder, New Zealand (45${}^{\circ }$ S, 170${}^{\circ }$ W). Measured number densities are shown as box plots (minimum excluding outliers below the 0.4th percentile, 25th percentile, median, 75th percentile, maximum excluding outliers above the 99.6th percentile) and modeled number densities as solid lines. For the Laramie plots (a), OPC measurements are used from the period 1999–2008 and model results are averaged over the 5 years of the time slice. For the Lauder plots (b), OPC measurements are used from January to April 1998–2001 and zonal mean model results are averaged from January to April over 5 years of the time slice. Model results are weighted with the counting efficiencies for OPC channels from for direct comparability with the measurements.

[Figure omitted. See PDF]

SOCOL-AERv2 simulates higher number densities of condensation nuclei (CN, $r>\mathrm{0.01}$ $\mathrm{\mu }\mathrm{m}$) above 25 km, matching the shape of the observed curve better than SOCOL-AERv1. The transport of polar ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$-rich air to midlatitudes during the breakup of the polar vortex may lead to the high number densities of CN above 25 km in the measurements . The improved agreement in SOCOL-AERv2 is due to the implementation of dry radius binning and the improvement in sulfur mass conservation, which enable the model to capture the increased transport of CN to midlatitudes during late winter and spring. SOCOL-AERv2 also displays a kink at the tropopause for particle channels larger than $r>\mathrm{0.15}$ $\mathrm{\mu }\mathrm{m}$, which appeared after the addition of interactive wet deposition. In the interactive wet deposition scheme, the scavenging efficiency of particles depends on radius, which leads to stronger removal of larger aerosol particles in the troposphere. Lauder OPC measurements may also show a similar kink at the tropopause for the larger particle channels; however, it is difficult to verify this given the large variability in the measurements (Fig. ).

Otherwise, both model versions show similar levels of agreement with the OPC measurements. All four channels larger than $r>\mathrm{0.25}$ $\mathrm{\mu }\mathrm{m}$ have too high number densities compared to observations at Laramie, with the agreement becoming even worse with altitude. At Lauder the agreement of the largest three channels ($r>\mathrm{0.30}$ $\mathrm{\mu }\mathrm{m}$) is better. attributed the worsening agreement of larger aerosol particles with altitude to either numerical diffusion in the sedimentation scheme or an overestimate in the speed of the Brewer-Dobson circulation, a known artifact in the SOCOL model .

3.2.4

Comparison with UTLS ${\mathrm{SO}}_{\mathrm{2}}$ measurements

There has been an ongoing debate in the literature regarding the magnitude of the cross-tropopause ${\mathrm{SO}}_{\mathrm{2}}$ flux and its relative importance in establishing the Junge layer . The debate has been fueled by a lack of in situ measurements in the UTLS region and the high temporal and spatial variability in UTLS ${\mathrm{SO}}_{\mathrm{2}}$. Figure  compares the tropical UTLS ${\mathrm{SO}}_{\mathrm{2}}$ measured by two aircraft campaigns with two annual mean satellite products, MIPAS , and ACE-FTS , averaged during volcanically quiescent periods. The two in situ measurements and the ACE-FTS satellite product all show ${\mathrm{SO}}_{\mathrm{2}}$ mixing ratios of 5 to 10 pptv around the tropical tropopause ($\sim \mathrm{17}$ km). MIPAS-observed ${\mathrm{SO}}_{\mathrm{2}}$ is around 24 pptv at 17 km, substantially higher than the other observations.

Figure 5

Comparing modeled and measured ${\mathrm{SO}}_{\mathrm{2}}$ mixing ratios in the tropical UTLS region between 10 and 25${}^{\circ }$ N. Modeled results from three simulations are averaged over 5 years of the time slice and are shown as colored dashed lines. Horizontal error bars indicate the interquartile range of monthly means. Observational datasets are shown as solid lines with interquartile ranges, extracted from data in Fig. 3 of and Fig. 11 of . The satellite datasets (triangles), MIPAS and ACE-FTS, show mean ${\mathrm{SO}}_{\mathrm{2}}$ values between 2002–2012 and 2004–2010, respectively, and have been filtered to remove any data affected by major volcanic eruptions. data (brown circles) represent in situ flight data from October 2015 over the Gulf of Mexico and the tropical eastern Pacific Ocean. data (orange circles) were measured in a flight campaign over the tropical western Pacific Ocean in October 2016. Model and satellite data are averaged over all longitudes between 10 and 25${}^{\circ }$ N.

[Figure omitted. See PDF]

The annual means of SOCOL-AERv1 and v2 agree with MIPAS-observed ${\mathrm{SO}}_{\mathrm{2}}$ and overestimate the three other observation sets at the tropopause, simulating ${\mathrm{SO}}_{\mathrm{2}}$ mixing ratios between 20 and 30 pptv at 17 km. MIPAS satellite observations of ${\mathrm{SO}}_{\mathrm{2}}$ under nonvolcanic conditions are uncertain, which may explain the systematic offset between MIPAS ${\mathrm{SO}}_{\mathrm{2}}$ measurements and the other observation sets . On the other hand, in situ measurements lack the spatial and temporal coverage of satellites, which reduces their comparability with global models. More aircraft campaigns will be invaluable for determining the background level of UTLS ${\mathrm{SO}}_{\mathrm{2}}$. If anything, the currently available observations suggest that SOCOL-AER's cross-tropopause ${\mathrm{SO}}_{\mathrm{2}}$ transport might be too high. In Sect.  we will investigate the consequences of an overestimated UTLS ${\mathrm{SO}}_{\mathrm{2}}$.

3.3 Observational disagreements with SOCOL-AERv2

To summarize, SOCOL-AERv2 shows similar levels of agreement with stratospheric sulfur observations as SOCOL-AERv1. The stratospheric aerosol burden simulated by SOCOL-AERv2, 160 Gg S, agrees very well with the SAGE-$\mathrm{3}\mathit{\lambda }$ retrieved burden, 165 Gg S. SOCOL-AERv2 slightly overestimates the SAGE II aerosol extinction in the lowermost stratosphere at the Equator and in the lowermost stratosphere at 45${}^{\circ }$ N, namely by up to 25 % at wavelength 1020 nm and by up to 40 % at 525 nm. Since organic particles may contribute to the aerosol burden in the lowest 1–3 km of the stratosphere, we think that SOCOL-AERv2 is actually overestimating the cross-tropopause transport of sulfur. OPC measurements also show that large particle channels ($r>\mathrm{0.25}$ $\mathrm{\mu }\mathrm{m}$) are overestimated in the UTLS at midlatitudes (by up to a factor of 3). A second region where SAGE II extinctions diverge from the simulated values is between 25 and 30 km at the Equator, where SOCOL-AERv2 overestimates extinctions. Above 35 km at the Equator and above 30 km at 45${}^{\circ }$ N, SOCOL-AERv2 underestimates aerosol extinctions; however, this is likely caused by a lack of meteoritic material in the model.

To address the possible overestimation of sulfur transport to the stratosphere, we ran two additional simulations, AER-SCAV and ICE-OX. In AER-SCAV, we increased the scavenging coefficient of aerosol on ice clouds by a factor of 20, from 0.05 to 1, maximizing the effect of upper tropospheric sulfate aerosol removal. In ICE-OX, the ice water content was added to the liquid water content before the aqueous-phase chemistry routine so that ${\mathrm{SO}}_{\mathrm{2}}$ oxidation occurs as well in middle and upper tropospheric ice clouds, maximizing the effect of condensed phase S(IV) to S(VI) oxidation. One laboratory study has identified the ${\mathrm{SO}}_{\mathrm{2}}$ $+$ ${\mathrm{H}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{2}}$ reaction on the surface of ice as a possible sink for ${\mathrm{SO}}_{\mathrm{2}}$, although complicating factors, like partial-pressure-dependent reaction probabilities and surface poisoning during the reaction, make it difficult to extrapolate the measurements to atmospheric conditions . Furthermore, found improved model agreement with Arctic sulfate measurements when they included ${\mathrm{SO}}_{\mathrm{2}}$ oxidation also in ice water. In addition, physical uptake of ${\mathrm{SO}}_{\mathrm{2}}$ without conversion to S(VI) on ice has been observed in the laboratory and may lead to gravitational settling; uptake of ${\mathrm{SO}}_{\mathrm{2}}$ on ice is not considered in either SOCOL-AERv1 or v2. Assuming that ${\mathrm{SO}}_{\mathrm{2}}$ oxidation occurs in cloud ice water at the same rate as cloud liquid water is likely an upper limit estimate for the scavenging of ${\mathrm{SO}}_{\mathrm{2}}$ on ice.

These extreme simulations succeed in reducing the cross-tropopause sulfur transport, leading to strongly reduced stratospheric aerosol burdens, namely 133 Gg S in AER-SCAV and 92 Gg S in ICE-OX. In these two simulations, extinctions at 45${}^{\circ }$ N now match observations in the lowermost stratosphere, while equatorial extinctions underestimate observations, which may be reasonable since organic aerosol particles play a role in this level (Fig. S2). Similarly, the modeled OPC channels are reduced in number density in AER-SCAV and ICE-OX at midlatitudes (Fig. S3). However, ICE-OX now shows too low number densities of CN in the UTLS compared to OPC measurements, suggesting either that too much ${\mathrm{SO}}_{\mathrm{2}}$ is removed or that other aerosol types contribute to CN at these altitudes. The available ${\mathrm{SO}}_{\mathrm{2}}$ measurements also imply that too much ${\mathrm{SO}}_{\mathrm{2}}$ is removed in ICE-OX since the simulated ${\mathrm{SO}}_{\mathrm{2}}$ concentration at 17 km ($\sim \mathrm{1}$ pptv) is lower than the in situ and ACE-FTS values of 5 to 10 pptv (Fig. ).

It is important to mention that although the agreement in the UTLS was improved by AER-SCAV and ICE-OX, there may be other reasons behind the too high cross-tropopause transport in SOCOL-AERv2. Convective transport of ${\mathrm{SO}}_{\mathrm{2}}$ and aerosol to the upper troposphere may be too strong in SOCOL-AERv2, which is a common problem with other GCMs . In this case, chemical oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ on ice or increased aerosol scavenging, rather than being a missing feature in itself, would be compensating for the strong convective transport. Numerical diffusion may further enhance the ${\mathrm{SO}}_{\mathrm{2}}$ cross-tropopause transport due to the strong vertical gradient at the tropopause. The influence of convection could be further investigated by testing the sensitivity of SOCOL-AER's sulfur cycle to different convection schemes, as has been done for other models . The choice of the convective scheme and the order in which it is called relative to the wet deposition routine could be used to further tune SOCOL-AERv2 in the UTLS; however, a clear challenge is the lack of in situ measurements at these altitudes. Further measurements of these species in the UTLS would be helpful to constrain the importance of ${\mathrm{SO}}_{\mathrm{2}}$ and primary sulfate aerosol in establishing the stratospheric aerosol layer.

3.4 Evaluation of SOCOL-AER deposition in transient simulations

In order to evaluate the performance of SOCOL-AER versions in the troposphere, we will compare simulated annual mean deposition fluxes with the database compiled for the World Meteorological Organization (WMO) assessment of precipitation chemistry and deposition . The WMO assessment only included regionally representative sites, e.g., excluding measurements within 50 km of industrial or urban areas, which should be comparable to the simulated values in a global model with coarse resolution. The deposition fluxes reported in the WMO assessment were averaged in 3-year periods, 2000–2002 and 2005–2007. The WMO assessment corrected wet deposition measurements for sea salt contributions of sulfate at sites less than 100 km from coastlines and at all African measurement sites. For our study, we only use the sites where measurement methodology and temporal data coverage were assessed as “satisfactory” or “conditional” in the WMO database. We interpolated modeled annual mean deposition to the coordinates of the measurements stations. Several previous model intercomparison projects that simulated deposition will be used as benchmarks for the performance of SOCOL-AER compared with observations.

3.4.1 Wet deposition

The wet deposition scheme in SOCOL-AERv2 is coupled to the climate model’s cloud and precipitation fields, whereas in SOCOL-AERv1 constant wet deposition lifetimes are applied. SOCOL-AERv1 therefore simulates too large deposition fluxes in dry regions and too small deposition fluxes in wet regions compared to SOCOL-AERv2 (Sect. ). To verify the SOCOL-AERv2 wet deposition fluxes in both dry and wet regions, we compare simulated and observed sulfate deposition fluxes over 3 orders of magnitude in Fig. . The presentation is logarithmic since using linear axes would give too high a weight to large deposition fluxes, obscuring biases at the lower range of deposition fluxes.

Figure 6

Evaluation of modeled total sulfur wet (a, b, e, f) and dry (c, d, g, h) deposition fluxes against measurement sites from the WMO database . SOCOL-AERv1 and SOCOL-AERv2 are compared with measurements averaged in two different time periods, 2000–2002 and 2005–2007. The ensemble standard deviation for the model results is shown as vertical bars. A power regression between the simulation results and measurements is shown in blue and can be compared to the one-to-one line shown in black. Two model evaluation metrics are listed on the plots: the goodness of fit of the power regression between model and measurements ($R^{\mathrm{2}}$) and the fraction of stations where the model is within a factor of 2 of measurements ($f_{\mathrm{2}\times }$). Points are colored according to the region (for the wet deposition plots) or the measurement network (for the dry deposition plots) of the measurement stations.

[Figure omitted. See PDF]

SOCOL-AERv1 indeed overestimates low deposition fluxes ($<\mathrm{1}$ kg S ha${}^{-\mathrm{1}}$), corresponding to sites in drier areas ($<\mathrm{50}$ cm yr${}^{-\mathrm{1}}$ precipitation). For both time periods (2000–2002 and 2005–2007), SOCOL-AERv2 improves the agreement with observations compared to SOCOL-AERv1 ($R^{\mathrm{2}}=\mathrm{0.61}$ and $R^{\mathrm{2}}=\mathrm{0.69}$ for SOCOL-AERv2 vs. $R^{\mathrm{2}}=\mathrm{0.51}$ and $R^{\mathrm{2}}=\mathrm{0.58}$ for SOCOL-AERv1). The fraction of stations where the model is within a factor of 2 of observations ($f_{\mathrm{2}\times }$) also improves slightly for both measurement periods in SOCOL-AERv2. The variability in simulated wet deposition fluxes, shown by the ensemble standard deviation bars in Fig. , increases in SOCOL-AERv2 because wet deposition is coupled to modeled precipitation. Due to the internal variability in modeled precipitation in a free-running climate model, multiple-ensemble simulations as well as long-term deposition measurements are required when comparing models with observations. Overall, SOCOL-AERv2 matches the measurements better than SOCOL-AERv1, especially for sites with low and high deposition fluxes.

Nevertheless, there are remaining biases in the deposition fields of SOCOL-AERv2, e.g. high biases in many North American sites compared to the WMO observations (Fig.  and Table ). Since the model's deposition scheme is coupled to precipitation fields, inaccuracies in the modeled precipitation distribution can lead to incorrect deposition fluxes. We calculated the model's precipitation biases compared to the WMO database for each station. The bias in precipitation depth in SOCOL-AERv2 correlates with the bias that we find for sulfate deposition fluxes (Spearman's correlation coefficient, $\mathit{\rho }=\mathrm{0.5}$, Fig. S4). This finding suggests that some of the model biases can be explained by errors in precipitation fields rather than errors in the wet deposition scheme. Both versions of SOCOL-AER match the observations better in the period 2005–2007 rather than 2000–2002. Since there are no large differences in the precipitation biases between these periods, this improvement could be related to more accurate ${\mathrm{SO}}_{\mathrm{2}}$ emission maps for the 2005–2007 period. One known issue with the CEDS anthropogenic ${\mathrm{SO}}_{\mathrm{2}}$ inventory is that the emissions from the western United States are overestimated compared to the eastern United States . SOCOL-AERv2 also shows higher deposition biases in the western United States. Since errors in emission inventories and model precipitation fields impact the evaluation of the deposition field, it is difficult to attribute errors to the deposition scheme itself.

Table 3

Metrics comparing the agreement of simulated sulfate wet deposition with the WMO database, separated by measurement region. Results from SOCOL-AERv2 from two time periods are compared with the range of values from past model intercomparison projects (MIP), including Photocomp , HTAP I , ACCMIP , and HTAP II . Note that the observational and simulation time periods covered by the other model intercomparison projects differ from SOCOL-AER in this study (see Table S1 in the Supplement). The metrics are calculated in linear space to conform with past MIPs: linear fit slopes differ from the power regression in Fig.  and fractions of sites where the model is within $\pm \mathrm{50}$ % of observations differ from the fraction within a factor of 2 listed in Fig. .

Metric	North America			Europe			East Asia
	Range of MIPs	SOCOL-AERv2		Range of MIPs	SOCOL-AERv2		Range of MIPs	SOCOL-AERv2
		2000 to	2005 to		2000 to	2005 to		2000 to	2005 to
		2002	2007		2002	2007		2002	2007
Mean observations	2.5–3.1	3.3	3.1	2.3–4.0	3.8	2.9	6.5–6.9	8.3	7.8
(kg S ha${}^{-\mathrm{1}}$)
Mean model	2.8–3.6	4.4	4.0	2.0–4.6	3.7	3.5	3.9–5.0	5.8	5.9
(kg S ha${}^{-\mathrm{1}}$)
Mean bias	$-$0.2 to 0.5	1.1	0.9	$-\mathrm{1.3}$ to 0.5	$-\mathrm{0.1}$	0.6	$-\mathrm{2.9}$ to $-\mathrm{1.6}$	$-\mathrm{2.5}$	$-\mathrm{1.9}$
(kg S ha${}^{-\mathrm{1}}$)
Linear fit slope	0.6–1.0	1.2	1.1	0.3–0.6	0.5	0.6	0.3–0.5	0.4	0.4
(observations vs.
model)
Pearson correlation	0.8–0.9	0.8	0.8	0.6–0.7	0.6	0.5	0.6–0.9	0.8	0.7
coefficient, $R$
Fraction of sites	70–77	48	60	53–86	72	61	69–88	66	73
within $\pm \mathrm{50}$ %
(%)

Table  compares the performance of SOCOL-AERv2 to past model intercomparison studies, including Photocomp , ACCMIP , HTAP I , and HTAP II . These intercomparison projects used observational data from the same networks as in the WMO database to evaluate their results. However, the analysis periods for these intercomparison projects differed from this study for both the simulations and observations (Table S1), which can contribute to the differences in the results. SOCOL-AERv2 shows similar levels of agreement with observations as the previous model intercomparison studies in the European and East Asian regions. The model biases, correlation coefficients, and fraction of values within $\pm \mathrm{50}$ % of measurements fall within or very close to the range of the intercomparison projects. In North America, SOCOL-AERv2 correlates similarly with observations (Pearson correlation coefficient, $R=\mathrm{0.8}$); however, biases and the fraction within $\pm \mathrm{50}$ % are worse than the past intercomparison projects. As mentioned above, the North American deposition fluxes may be affected by inaccuracies in our model's precipitation fields and/or errors in the anthropogenic ${\mathrm{SO}}_{\mathrm{2}}$ emission inventory. In addition, there are several factors that advantage the model-intercomparison projects compared to the SOCOL-AER runs. Firstly, we have compared ranges of the multimodel means and not the individual model values from the intercomparison projects, which likely have a wider spread and individually worse performance. Secondly, model resolution can play an important role in the comparison between observations and simulations and SOCOL-AER was run at a relatively coarse resolution ($\mathrm{2.8}{}^{\circ }\times \mathrm{2.8}{}^{\circ }$) compared to many of the models used in the intercomparison project. Finally, three of the intercomparison projects (Photocomp, HTAP I, and HTAP II) were based on offline chemistry-transport models that are run with observed meteorology. On the other hand, SOCOL-AER is used in free-running mode, producing five ensemble members that are subsequently averaged. SOCOL-AER, and its parent climate model ECHAM5, include precipitation biases that impact the simulation of deposition. Considering these aspects, the performance of SOCOL-AERv2 compares well with state-of-the-art sulfate wet deposition models.

3.4.2 Dry deposition

Dry deposition fluxes compiled by the WMO assessment are based on North American stations from the Canadian Air and Precipitation Monitoring Network (CAPMoN) and the US Clean Air Status and Trends Network (CASTNET). Dry deposition fluxes are not measured directly but are inferred through surface-based measurements of gas and particle concentrations and estimates of their dry deposition velocities . The estimation of dry deposition velocities is uncertain; the ${\mathrm{SO}}_{\mathrm{2}}$ dry deposition velocities calculated by the CAPMoN network are around 50 % higher than those of the CASTNET network . The inferred dry deposition fluxes are therefore less reliable for comparisons with models than the measured wet deposition fluxes.

Figure  displays the comparison between SOCOL-AER versions and the total sulfur (sum of ${\mathrm{SO}}_{\mathrm{2}}$ and aerosol) inferred dry deposition fluxes from North American measurement networks. Despite the systematic bias between the CAPMoN and CASTNET networks, we calculate the power regression, $f_{\mathrm{2}\times }$, and the $R^{\mathrm{2}}$ value based on the combined dataset of both networks. We take this approach as it is unclear which network's dry deposition velocity calculation is more accurate . SOCOL-AERv2, with its new interactive dry deposition scheme, shows a similar correlation with the observations as SOCOL-AERv1. However, SOCOL-AERv2's improved deposition scheme simulates much more realistic slopes of the correlation lines and higher fractions of model results within a factor of 2 of observations. The improved agreement of SOCOL-AERv2 is likely caused by the better spatial representation of ${\mathrm{SO}}_{\mathrm{2}}$ dry deposition velocities in the interactive scheme. The modeled dry deposition fluxes can also be affected by the new wet deposition scheme since enhanced sulfur removal through wet deposition leaves less sulfur available for dry deposition, and vice versa. Similar to the wet deposition comparison, the model performs better in 2005–2007 compared to 2000–2002, suggesting that the North American emission inventories may be biased in the earlier time period. Additional comparisons for the ${\mathrm{SO}}_{\mathrm{2}}$ and aerosol dry deposition fluxes with observations show better agreement for SOCOL-AERv2 than SOCOL-AERv1 (Fig. S5 and S6).

Compared to past model-intercomparison projects for total sulfur dry deposition, SOCOL-AERv2 simulates a lower bias and a higher fraction of model values within $\pm \mathrm{50}$ % of the observation sites (Table ). The agreement nevertheless remains poor with only 19 % to 28 % of the modeled total dry deposition fluxes within $\pm \mathrm{50}$ % of the observations. Aerosol dry deposition biases are larger in SOCOL-AER compared to past model intercomparison projects. However, since ${\mathrm{SO}}_{\mathrm{2}}$ dominates the dry deposition flux (compare Figs. S5 and S6), the reduced biases in the ${\mathrm{SO}}_{\mathrm{2}}$ deposition fluxes in SOCOL-AERv2 leads to overall lower total sulfur dry deposition biases. It is difficult to conclude whether the observations should actually be the target for the model in this case since dry deposition fluxes are inferred values with inherent uncertainties. As observational networks improve their parametrizations for deriving the dry deposition flux, they will become more reliable standards with which to compare modeled results. However, the similar, if not better, agreement of SOCOL-AERv2 compared to sulfur dry deposition from model-intercomparison studies adds confidence to the implemented dry deposition scheme.

Table 4

Metrics comparing the agreement of simulated dry deposition with the WMO database in North America, separated by measurement quantity (total sulfur, ${\mathrm{SO}}_{\mathrm{2}}$, and sulfate aerosol). Results from SOCOL-AERv2 from two time periods are compared with the range of values from past model intercomparison projects (MIP), including HTAP I , ACCMIP , and HTAP II . Note that the observational and simulation time periods covered by the other model intercomparison projects differ from SOCOL-AER in this study (see Table S1). The metrics are calculated in linear space to conform with past MIPs: linear fit slopes differ from the power regression in Fig.  and fractions of sites where the model is within $\pm \mathrm{50}$ % of observations differ from the fraction within a factor of 2 of observations listed in Fig. .

Metric	Total sulfur dry deposition			${\mathrm{SO}}_{\mathrm{2}}$ dry deposition			Sulfate aerosol dry deposition
	Range of MIPs	SOCOL-AERv2		Range of MIPs	SOCOL-AERv2		Range of MIPs	SOCOL-AERv2
		2000 to	2005 to		2000 to	2005 to		2000 to	2005 to
		2002	2007		2002	2007		2002	2007
Mean observations	1.1–2.3	2.5	2.0	0.8–1.9	2.0	1.6	0.2–0.4	0.4	0.4
(kg S ha${}^{-\mathrm{1}}$)
Mean model	3.6–5.1	4.9	3.8	3.2–4.6	3.8	2.9	0.4–0.5	1.1	0.9
(kg S ha${}^{-\mathrm{1}}$)
Mean bias	2.5–3.7	2.4	1.7	2.4–2.6	1.7	1.2	0.1–0.2	0.7	0.5
(kg S ha${}^{-\mathrm{1}}$)
Linear fit slope	1.0–2.7	1.3	1.1	1.0–2.7	1.1	0.9	1.0–1.6	2.1	2.0
(observations vs.
model)
Pearson correlation	0.8	0.8	0.8	0.8	0.8	0.8	0.8–0.9	0.6	0.6
coefficient, $R$
Fraction of sites	6–13	19	28	5–6	22	31	47–48	19	29
within $\pm \mathrm{50}$ %
(%)

3.5 Pinatubo simulation with SOCOL-AERv2

The eruption of Mt. Pinatubo in 1991 remains the strongest directly observed volcanic event, which makes it a valuable test for models. analyzed the performance of SOCOL-AERv1 for this case in a series of sensitivity experiments, demonstrating reasonable agreement with observations of different aerosol parameters. Figure  shows the SOCOL-AERv1 and v2 global total stratospheric aerosol burden for emission estimates of 6 and 7 Tg S compared to HIRS, SAGE-$\mathrm{3}\mathit{\lambda }$ (v4), and SAGE-$\mathrm{4}\mathit{\lambda }$. The same Mt. Pinatubo emission estimates modeled with SOCOL-AERv2 show clear differences compared to SOCOL-AERv1, expressed in both the shape of the peak and its magnitude. SOCOL-AERv2 is $\sim \mathrm{0.4}$ Tg S higher at the peak values in late 1991 but lower at the tail after mid-1992. The main reasons for a narrower and stronger peak in SOCOL-AERv2 are the changes to the chemistry scheme for the CCMI simulation and the update of reaction coefficients to recommendations (Sect. ), which led to higher OH concentrations in the UTLS. In SOCOL-AERv2, the oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ is therefore faster, causing faster aerosol formation and its earlier removal from the stratosphere. Another change that contributed to the differences is the improved sulfur mass conservation. This effect is illustrated by the integrated total sulfur deposition in the two model versions, after the emission of 7 Tg S with all other sulfur emissions switched off (orange and grey lines in Fig. ). This experiment shows that by the end of 1995, when almost all the volcanic aerosol is already removed from the atmosphere, SOCOL-AERv1 deposited only 6.4 Tg S, while in SOCOL-AERv2 this is now improved to 6.7 Tg S, i.e., the volcanic sulfur mass loss decreased from 8.6 % to 4.2 %. The remaining part of the mass loss is likely due to limitations of the transport scheme . Overall, in SOCOL-AERv2 the Mt. Pinatubo emission estimate required for the model to reproduce the observed burden peak is still between 6 and 7 Tg S, given the large observational uncertainty.

Figure 7

Time evolution of the globally averaged stratospheric aerosol burden calculated by SOCOL-AERv1 (Fig. 1a from ) and SOCOL-AERv2 compared with the HIRS and SAGE II-derived data (SAGE-$\mathrm{4}\mathit{\lambda }$ and SAGE-$\mathrm{3}\mathit{\lambda }$v4). Light blue shaded area: uncertainties in HIRS. Black shaded area: $\mathrm{2}\mathit{\sigma }$ five-member ensemble spread of one of the model experiments; others are shown as ensemble means. Model experiments are performed with two emission estimates: 6 and 7 Tg S. Orange and grey lines represent the accumulated (acc) globally integrated deposition of sulfur emitted from Pinatubo. The accumulated deposition from Pinatubo is calculated from simulations where all sulfur sources other than Pinatubo are turned off.

[Figure omitted. See PDF]

3.6 Updated nonvolcanic sulfur budget for the year 2000

Figure  shows an update from SOCOL-AERv2 for the atmospheric sulfur budget under volcanically quiescent background conditions. The model was run with repeating boundary conditions for the year 2000, with the updated sulfur emission datasets (Sect. ). Due to higher emissions in the updated sulfur emission datasets, the stratospheric aerosol burden shown in Fig.  (167 Gg S) is slightly larger than the burden for the SOCOL-AERv2 run in Table  (160 Gg S). We have added diagnostics to track tracer fluxes within the model's planetary boundary layer (PBL), expanding from the budget figure from SOCOL-AERv1 (Fig. 3 in ). For these calculations, we extracted the modeled height of the PBL from the ECHAM5 vertical diffusion routines. In Fig. , italicized fluxes are calculated from the other outputted fluxes assuming steady-state conditions, i.e., that the fluxes into/from a species add up to zero. The fluxes in the PBL do not fully balance due to failures in the steady-state assumption, difficulties in extracting chemical fluxes from the iterative chemical solver, or remaining mass conservation errors in the model. However, the imbalances in the PBL are at most 4 % of the total input fluxes in the PBL for each species. Therefore, the outputted fluxes shown in Fig.  can be considered reliable.

Figure 8

Atmospheric sulfur budget from SOCOL-AERv2 under volcanically quiescent conditions for the year 2000. The figure is produced using output diagnostics that track sulfur fluxes and burdens within the planetary boundary layer (PBL), free troposphere, and stratosphere. Solid arrows show net fluxes and dashed arrows show one-way fluxes, all in gigagrams of sulfur per year (Gg S yr${}^{-\mathrm{1}}$). Simulated burdens of sulfur species are given within the boxes, in units of gigagrams of sulfur (Gg S). Italicized numbers represent fluxes that are derived by balancing other fluxes, assuming steady state of the upper layers. For example, net cross-tropopause fluxes are calculated by balancing the stratospheric chemical fluxes and net cross-PBL fluxes are calculated by balancing the free tropospheric chemical fluxes. Upward OCS and sulfate aerosol cross-tropopause fluxes are calculated based on the residual meridional and vertical air velocities ($v^{*}$ and $w^{*}$) and concentrations at the tropopause. Black numbers: SOCOL-AERv2 results. Red numbers: (a) SAGE-3$\mathit{\lambda }$ stratospheric aerosol burden (b) DMS emissions estimated by (c) dry deposition fluxes from NCAR-CAM3.5 model for the year 2000, which participated in ACCMIP (d) multimodel mean wet deposition from ACCMIP models for the year 2000 .

[Figure omitted. See PDF]

Figure  reveals differences between the tropospheric sulfur cycle in the PBL and the free troposphere. OCS, MSA, and sulfate aerosol all show higher burdens in the free troposphere compared to the PBL since these species are produced chemically in the atmosphere and/or they are long-lived species. On the other hand, several sulfur species (CS${}_{\mathrm{2}}$, DMS, ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$, and ${\mathrm{SO}}_{\mathrm{2}}$) have higher burdens in the boundary layer compared to the free troposphere, despite the free troposphere having more volume. Emission from the surface is the most important atmospheric source for these species. Still, free tropospheric ${\mathrm{SO}}_{\mathrm{2}}$ is supplied not only by cross-PBL transport (67 %) but also through transport and oxidation of short-lived species (mainly DMS and ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$). This result supports another modeling study that highlighted the transport pathway of DMS to the upper troposphere by deep convection over the ocean . Aqueous-phase ${\mathrm{SO}}_{\mathrm{2}}$ oxidation dominates oxidation of ${\mathrm{SO}}_{\mathrm{2}}$ in both parts of the troposphere in SOCOL-AERv2 (74 % of total ${\mathrm{SO}}_{\mathrm{2}}$ oxidation in the PBL and 68 % in the free troposphere). Sulfate aerosol is produced in the PBL through condensation of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ on existing particles and aqueous-phase oxidation of ${\mathrm{SO}}_{\mathrm{2}}$, with negligible nucleation fluxes. In the upper troposphere, nucleation of new ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$-${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ droplets becomes more important, although growth of existing particles remains the largest mass flux to the particle phase. Figure  shows a large downward flux of aerosol from the free troposphere to the PBL (18 Tg S yr${}^{-\mathrm{1}}$). This balanced flux does not represent only gravitational sedimentation to the PBL but rather mostly wet removal of free tropospheric aerosol to the surface. The PBL diagnostics in SOCOL-AERv2 are a useful tool for understanding transport and transformation of sulfur species in the troposphere.

SOCOL-AERv1 calculated a total flux to stratospheric aerosol of 181 Gg S yr${}^{-\mathrm{1}}$ by summing the net cross-tropopause fluxes of gaseous sulfur species and the upward cross-tropopause flux of primary sulfate aerosol . With the modifications that were made in this paper, SOCOL-AERv2 now simulates a flux of 228 Gg S yr${}^{-\mathrm{1}}$ into stratospheric aerosol. Of this flux, 46 % is due to upward transport of primary tropospheric aerosol, 27 % due to ${\mathrm{SO}}_{\mathrm{2}}$, 21 % due to OCS, 4 % due to DMS, and 1 % due to ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$. These contributions are very similar to the contributions reported for SOCOL-AERv1 . As discussed in Sect. , the AER-SCAV and ICE-OX runs suggest that too much ${\mathrm{SO}}_{\mathrm{2}}$ and tropospheric aerosol are transported across the tropopause in SOCOL-AERv2. Future improvements of convection schemes and increased availability of reliable observational data will further constrain the accuracy of the ${\mathrm{SO}}_{\mathrm{2}}$ contribution to the stratospheric aerosol layer.

4 Conclusions

For SOCOL-AER to be used to study the tropospheric sulfur cycle, as well as the deposition response to volcanic eruptions, we implemented new features and applied several corrections to the code. Compared to SOCOL-AERv1, the implemented interactive deposition schemes in SOCOL-AERv2 result in much improved agreement with measurements from sulfur deposition networks. Our tests for Pinatubo showed that SOCOL-AERv2 now gives more aerosol mass in 1991 due to faster ${\mathrm{SO}}_{\mathrm{2}}$ oxidation. Better sulfur mass conservation allowed us to decrease the sulfur mass loss after Pinatubo from 8.6 % to 4.2 %. With the improved mass conservation in SOCOL-AERv2, we are also able to separate free tropospheric fluxes from PBL fluxes in the atmospheric sulfur budget, revealing that short-lived sulfur species (DMS and ${\mathrm{H}}_{\mathrm{2}}\mathrm{S}$) contribute strongly to ${\mathrm{SO}}_{\mathrm{2}}$ in the free troposphere. With respect to stratospheric aerosol observations, SOCOL-AERv2 shows similar levels of agreement as SOCOL-AERv1. The modeled estimate for the burden of background stratospheric aerosol has increased from 116 Gg S in SOCOL-AERv1 to 160 Gg S in SOCOL-AERv2. At the same time, the burden derived from SAGE extinctions and OPC measurements increased from the SAGE-$\mathrm{4}\mathit{\lambda }$ estimate of 117 Gg S by 40 % to the most recent estimate from the SAGE-$\mathrm{3}\mathit{\lambda }$ dataset of 165 Gg S. Given the uncertainty in the SAGE-$\mathrm{3}\mathit{\lambda }$ and SAGE-$\mathrm{4}\mathit{\lambda }$ estimates, this might be to some degree fortuitous. Therefore, it is more reliable to compare the model with the satellite extinction measurements directly. Aerosol extinctions in the lower stratosphere are overestimated by SOCOL-AERv2 more than SOCOL-AERv1. We speculate that this might be related to inaccuracies in the model's convective transport scheme, which were also present in SOCOL-AERv1 but compensated by double counting of the ${\mathrm{SO}}_{\mathrm{2}}$ aqueous oxidation flux. Nevertheless, uncertain processes related to ${\mathrm{SO}}_{\mathrm{2}}$ and aerosol scavenging by ice clouds could also lead to overestimation of lower stratospheric sulfate aerosol. Comparison with other models in the background experiments of the Interactive Stratospheric Aerosol Model Intercomparison Project (ISA-MIP) may help to identify the cause of SOCOL-AERv2's overestimated lower stratospheric aerosol.

The model developments presented here increase the applicability of SOCOL-AER to scientific questions in both the troposphere and stratosphere. Namely, due to its improved deposition fluxes, SOCOL-AERv2 is more suitable for comparison with ice-core-derived magnitudes of past volcanic eruptions and their atmospheric impacts e.g.,; modeling the atmospheric budget of cosmogenic isotopes, which attach to sulfate aerosols ; and studying future changes to sulfur deposition relevant to agriculture and ecosystems e.g.,. With its updated chemistry and improved sulfur mass conservation compared to SOCOL-AERv1, SOCOL-AERv2 is more reliable for studying the impacts of volcanic eruptions and stratospheric sulfate geoengineering.

Code and data availability

Since the SOCOL-AER code is based on ECHAM5, users must first sign the ECHAM5 license agreement before accessing the SOCOL-AER code (http://www.mpimet.mpg.de/en/science/models/license/, last access: 22 August 2019). SOCOL-AER code is then freely available upon request from the authors. The simulation data presented in this paper are available at 10.3929/ethz-b-000342078 . SAGE II data from the GloSSAC database can be found online at 10.5067/GloSSAC-L3-V1.0. The OPC measurements from the University of Wyoming were obtained at ftp://cat.uwyo.edu/pub/permanent/balloon/Aerosol_InSitu_Meas/ (last access: 22 August 2019). Deposition flux measurements can be downloaded online from the World Data Centre for Precipitation Chemistry at http://www.wdcpc.org/global-assessment-data (last access: 22 August 2019).

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-12-3863-2019-supplement.

Author contributions

AF conducted most of the simulations, analyzed the data, and drafted the paper. TS conducted and analyzed the simulations for the post-Pinatubo period. AF was supervised directly by LW, TP, and AS during the model development work. All authors participated in the model development, discussions about the results, and revisions to the article.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

Thanks to Jianxiong Sheng for providing plotting codes to analyze SOCOL-AER results and for providing further information about SOCOL-AERv1. Thanks to Debra Weisenstein for discussions about the original AER code. Thanks to Jean-François Lamarque for providing the sulfur deposition data for the NCAR-CAM3.5 model. We acknowledge the researchers behind the GloSSAC SAGE II data and the OPC measurements for use of their data. We also gratefully acknowledge the sources of precipitation chemistry and deposition data acknowledged on p. 92 of , Atmos. Environ., 93, 10.1016/j.atmosenv.2013.10.060. Color tables from ColorBrewer 2.0 (http://colorbrewer2.org/#t) were used for the figures in this paper.

Financial support

This research has been supported by the ETH Zurich (grant no. ETH-39 15-2). Timofei Sukhodolov and Eugene Rozanov have been supported from the Swiss National Science Foundation (VEC (grant no. 200021-169241) and POLE (grant no. 200020_182239)).

Review statement

This paper was edited by Samuel Remy and reviewed by Claudia Timmreck and one anonymous referee.