2.1.1 Ocean model MPIOM

The oceanic part, MPIOM, is formulated on an Arakawa-C grid in the horizontal and on z levels in the vertical direction and solves the primitive equations with the hydrostatic and Boussinesq approximations (Jungclaus et al., 2006,, 2013). Subgrid-scale parameterizations include lateral mixing on isopycnals and tracer transports by unresolved eddies. Vertical mixing is represented as a combination of the Richardson-number-dependent scheme and the wind-driven turbulent mixing in the mixed layer (for details, see Jungclaus et al., 2013). The horizontal grid is consistent with the MPI-ESM LR configuration, which implies a bipolar grid (GR1.5) featuring one grid pole under Greenland and one under Antarctica. The resolutions are then regionally enhanced in the deep water formation regions and the overflows across the Greenland–Scotland Ridge so that the grid varies between 22 and 350 km. In the vertical, 40 levels are unevenly placed in the water column, with the first 20 levels distributed over the top 700 m. The bottom topography is represented by a partial-step formulation (Wolff et al., 1997).

The sea ice model combines the codes of MPIOM and ECHAM. In ECHAM, a simplified thermodynamic sea ice model is incorporated to provide at each atmospheric time step a physically consistent surface temperature in ice-covered regions. This part also contains a melt-pond scheme, which divides the surface of the sea ice into snow, bare ice, and melt pond with individual albedos (Pedersen et al., 2009). The atmospheric part of the sea ice model then integrates all surface fluxes into ice and provides this information to the oceanic part of the code, which uses it to calculate the sea ice surface energy balance and related changes in ice thickness. The calculations of sea ice concentration and thickness are based on the Semtner (1976) formulation and tuned in MPI-ESM to produce the annual average pre-industrial Arctic sea ice volume of roughly 20 000–25 000 km³ (see Mauritsen et al., 2019). Sea ice dynamics is calculated following a viscoplastic approach of Hibler (1979).

2.1.2 Marine biogeochemistry model HAMMOC

Ocean biogeochemistry in MPI-ESM is represented by the Hamburg Ocean Carbon Cycle (HAMOCC) model (Ilyina et al., 2013; Paulsen et al., 2017). It simulates the oceanic cycles of carbon and other biogeochemical elements such as nutrients (phosphate, nitrate, and iron), oxygen, silicate, phytoplankton, zooplankton, and detritus. HAMOCC includes biogeochemical processes in the water column, the sediment, and at the air–sea interface. Biogeochemical tracers in the water column are fully advected, mixed, and diffused by the flow field of MPIOM. In total, the model has 17 state variables calculated prognostically in the water column and 12 state variables in the sediment. Nitrogen-fixing cyanobacteria was added to the model as an additional prognostic phytoplankton class by Paulsen et al. (2017).

2.1.3 Land surface model JSBACH

The Jena Scheme for Biosphere-Atmosphere Coupling in Hamburg (JSBACH) is the land component of MPI-ESM1.2. It provides the lower boundary conditions for the atmosphere over land and describes the dynamics of the land biogeochemistry in interaction with global climate. JSBACH treats processes like soil hydrology (five-layer scheme of Hagemann and Stacke, 2015), soil and litter decomposition, land use change (tiling approach with 12 plant functional types and two types of bare surface), fires, and a nitrogen cycle (Goll et al., 2017). Note that soil and marine chemical schemes are not yet combined with the atmospheric chemistry scheme in the current model version.

2.1.4 ECHAM6

ECHAM6 is an atmospheric general circulation model (GCM) that describes the large-scale circulation and its coupling to diabatic processes, both of which are ultimately driven by radiative forcing. It consists of a dry spectral-transform dynamical core, a transport model, and a suite of physical parameterizations for the representation of diabatic processes. The prognostic variables are temperature, vorticity, divergence, logarithm of surface pressure, and humidity, as well as cloud ice and water. Tracer transport and diabatic processes (also referred to as “model physics”) are calculated on a Gaussian transform grid. The adiabatic core of ECHAM6 consists of a mixed finite-difference/spectral discretization of the primitive equation that is identical to that employed in ECHAM5 (Stevens et al., 2013). All major changes relative to ECHAM5 are therefore related to the model physics. These changes include an improved representation of radiative transfer in the shortwave (or solar) part of the spectrum; a completely new description of tropospheric aerosol; an improved representation of surface albedo, including the treatment of melt ponds on sea ice (see Sect. 2.1.1); and an improved representation of the middle atmosphere through the gravity wave forcing. In addition, minor changes have been made in the representation of convective processes. Several coding errors in model physics were also corrected in the latest version (see details in Mauritsen et al., 2019).

Transport of species is performed with the flux-form semi-Lagrangian scheme of Lin and Rood (1996). Though this scheme is mass conservative by design, its application on the sigma–pressure coordinate system can cause a violation of the mass conservation especially in the case of large spatial gradients (Jöckel et al., 2001; Stenke et al., 2013). Turbulent mixing adopts an eddy diffusivity and viscosity approach following Brinkop and Roeckner (1995). Moist convection is parameterized according to Tiedtke (1989), with extensions by Nordeng (1994) and Möbis and Stevens (2012). Stratiform clouds are computed diagnostically based on a relative humidity threshold (Sundqvist et al., 1989). As in ECHAM5, gravity wave drag (GWD) is calculated using a subgrid orography scheme (Lott, 1999). The propagation and dissipation of the waves follow the formulation of Palmer et al. (1986) and Miller et al. (1989). Non-orographic GWD parameterizations are based on a wave-spectrum approach (Hines, 1997a, b). However, some parameters of both orographic and non-orographic gravity wave schemes have been adjusted for use in ECHAM6 during the tuning process of MPI-ESM (see Mauritsen et al., 2019). Both longwave and shortwave radiative-transfer calculations are now described by the PSrad scheme (Pincus and Stevens, 2013), which is based on the k-correlated method of the Rapid Radiative Transfer Model for GCMs (RRTM-G) (Iacono et al., 2008). Additional extra-heating parameterization in the Hartley, Huggins, and Schumann–Runge bands and the Lyman-α line is now applied for better representation of the solar cycle in the mesosphere and stratosphere (Sukhodolov et al., 2014). The optical properties for radiation are updated every 2 h. In contrast to the base model version, which applies climatological fields for this purpose, the radiation calculation of SOCOLv4 uses the prognostic tracer concentrations of sulfate aerosol and all radiatively active species (e.g., ozone and methane) except CO₂. Concentrations for CH₄, N₂O, and CFCs are prescribed only at the lowermost model level, while CO₂ is prescribed for the entire atmosphere. Note that for photolysis rates calculations, we use a separate subroutine based on the application of lookup tables (Rozanov et al., 1999). Cloud scattering is parameterized according to Mie theory using maximum-random cloud overlap and an inhomogeneity parameter to account for three-dimensional effects. Surface albedo is parameterized according to Brovkin et al. (2013). The radiative properties of tropospheric aerosols are now represented by the Max Planck Institute Aerosol Climatology version 2 simple plume implementation (MACv2-SP) parameterization (Stevens et al., 2017; Fiedler et al., 2017). MACv2-SP is formulated in terms of nine spatial plumes associated with different major anthropogenic source regions. It prescribes the aerosol optical depth, asymmetry factor and single scattering albedo as a function of geographical position, height above ground level, time, and wavelength. The evolution and distribution of anthropogenic aerosol is approximated with mathematical functions, while the natural aerosol is prescribed climatologically based on Kinne et al. (2013). MACv2-SP also prescribes aerosol–cloud interactions in the form of a Twomey effect (Twomey, 1977), thus allowing the increase in the cloud droplet number concentration and associated reduction in droplet size under constant liquid water path.

3.1.1 Reference data for evaluation

We validate the simulated climate variables from SOCOLv4 against observations. For this purpose, we use reanalysis data because they incorporate various observations combined in one dataset using comprehensive techniques that allow gaps in space and time to be estimated. It should be noted that the combination of different data can lead to some errors in the reanalysis data, both in trends and variability, that have a purely methodological nature, like those due to changes in the observation systems or aging of instruments as well as the details of underlying models used for assimilation. Intercomparison of different reanalysis products by Long et al. (2017) showed that reanalyses agree in the lower–middle stratosphere but deviate more from each other in the upper stratosphere and lower mesosphere because of the decreased availability of direct observations and thus the increased dependence on the assimilating model details, like the model top and vertical resolution. This should be considered when making conclusions about deviations between model results and reanalysis. To validate the thermodynamical structure of the modeled atmosphere, we use the recent ERA5.1 reanalysis dataset (Simmons et al., 2020; Hersbach et al. 2020) from ECMWF (European Centre for Medium-Range Weather Forecasts) and MERRA-2 (Modern-Era Retrospective Analysis for Research and Applications version 2, Gelaro et al., 2017). For surface air temperature, we additionally use the BEST data (the Berkeley Earth Surface Temperatures, Cowtan et al., 2015), which are a merged land–ocean dataset.

3.1.2 Temperature and winds

To assess the ability of SOCOLv4 to reproduce the climate and its response to historical forcing, we show variability in annual mean surface air temperature and its anomaly (Fig. 3a, b), the Arctic sea ice extension anomaly (Fig. 3c), as well as the El Niño–Southern Oscillation (ENSO) Niño3.4 statistics (Fig. 4a, b) for the historical period 1980 to 2018. Anomalies were calculated as a deviation from the mean over the 1980–2018 period. Figure 5a, c also show the global mean 2 m temperature and its anomaly simulated with the SOCOLv4 and various reanalysis data. The model reveals a warming trend similar to observations, whereas the ensemble mean curve shows less interannual variability, as expected (Fig. 3a). As can be seen from the figure, the reanalysis data slightly disagree in terms of the absolute level and the model data are within this spread. Modeled global mean surface air temperature is somewhat warmer than the ERA5.1 and MERRA-2 data but colder than the BEST data. Temperature anomalies (Fig. 3b) agree very well with observational composites. The annual Arctic sea ice extent anomalies are also accurately reproduced by SOCOLv4 in line with observations, showing a pronounced decline with time because of global warming and fluctuations of magnitudes to observations (Fig. 3c). According to Ding et al. (2019), the spread in sea ice extent between the ensemble members can be attributed to internal variability driven by fluctuations in Arctic pressure expressed as high pressure enhancing the sea ice loss and low pressure restraining it. Based on the analysis of Fig. 3, we conclude that the model response to historical radiative forcing is sufficiently close to observations, which proves suitability of the model for studies of past and future climate changes.

Figure 3Global and annual mean 2 m temperature (a) and its anomalies relative to 1980–2018 mean (b) and Arctic sea ice extent annual anomaly (c) in comparison to reanalysis data. The thick black line represents simulated mean of the three-member SOCOLv4 ensemble, thin black lines are individual ensembles members, the red line indicates ERA5.1 reanalysis, the green line indicates MERRA-2 reanalysis, and the blue line indicates BEST reanalysis.

Download

ENSO is known as the largest mode of sea surface temperature anomalies with irregular period of ∼ 2–7 years, which influences atmospheric circulation globally in the troposphere (e.g., Bulgin et al., 2020) as well as in the middle atmosphere (Domeisen et al., 2019). It is tracked as a variability of the sea surface temperature in the so-called Niño3.4 region (5^∘ N–5^∘ S, 120–170^∘ W). These temperature anomalies are usually in the range of ±3 K with positive fluctuations being associated with El Niño events and negative with La Niña events. There were several El Niño events observed with positive deviations stronger than +1.5 K during the considered decades, with the strongest being those in 1982/1983, 1997/1998, and 2015/2016. To verify the model performance in this context, we analyzed the fast Fourier transform (FFT) and the probability density function (PDF) (Fig. 4a and b, respectively) of the modeled and ERA5.1 Niño3.4 indices. The analysis indicates that SOCOLv4 has realistic ENSO in terms of amplitude and periodicity; however, the 40-year period is rather short for a robust estimation of ENSO frequencies, which is illustrated as a pronounced variability among the ensemble members. The common bias shared by all members is that ENSO events in the model mostly last longer than in observations. This is expressed as the power spectrum shift to shorter frequencies and a flattened PDF. It must be noted that this is a known feature of our underlying model, MPI-ESM, and our FFT results are very similar to those shown by Müller et al. (2018) in their Fig. 13, where they analyzed a 140-year period. The spatial pattern of ENSO in the model is reasonable but also has some biases related to the asymmetry between El Niño and La Niña events (for details see Tian et al., 2019). This drawback could be minimized by increasing coupling frequency between ocean and atmosphere from daily to hourly. However, this change also increases the computational demands to levels, which are less suitable for long-term climate calculations. Therefore, the LR version of MPI-ESM which we used for SOCOLv4 was tuned with a daily coupling (Mauritsen et al., 2019).

Figure 4(a) The 1980–2018 Niño3.4 power spectra (^∘C/cycles month⁻¹) calculated from the ERA5 data (orange) as well as from the three SOCOLv4 ensemble members (black). Dashed lines are the 95 % significance levels based on an auto-regressive AR(1) process fitting. (b) Histogram of Niño3.4 temperature anomalies.

Download

The left panel of Fig. 5 shows the geographical distribution of the 2 m air temperature climatology for 1980–2018 as simulated with SOCOLv4 and derived from the ERA5.1 reanalysis, as well as their difference. ERA5.1 data have been interpolated to the model grid. SOCOLv4 properly reproduces the observed geographical pattern of the surface air temperature with cold high-latitude areas in both hemispheres and warm tropics and subtropics. However, there are some regional discrepancies between the model and the reanalysis data. The model surface air temperature is generally slightly (∼ 0.2 K) warmer than in ERA5.1, which is also seen from the global means (Fig. 3a), but there are also several spots with cold biases. In the tropical area, the discrepancies are within ±1 K, except in the eastern part of South America over the Brazilian Highlands, where negative bias reaches −5 K, and over the Tibetan Plateau and mainland southeast Asia, with up to −3 K bias. In the northern middle latitudes, we obtain negative discrepancies of about −5 K over the Rocky Mountains and the North Atlantic Ocean. The most pronounced biases of ±6 K are located in the Southern Hemisphere over Antarctica. Most of these biases in high-altitude regions are a common feature of all coupled models, which is usually attributed to the simplified model's topography and uncertainties in the observational data (e.g., Guo et al., 2020). The right panel of Fig. 5 shows the 2 m surface air temperature trends over 1980–2018 over the Northern Hemisphere simulated by SOCOLv4 and calculated from ERA5.1 reanalysis data. The linear trend and its significance were calculated using a robust nonparametric Sen–Mann–Kendall test applying 95 % confidence interval. Simulated and observed trends over the Northern Hemisphere show a good spatial agreement including the location of the warming maximum north of Siberia. The maximum simulated trend amounts to 1.69±0.16 K, which is only about 0.25 K smaller than the maximum trend of 1.93 K derived from the reanalysis.

Figure 5Left panels: global distribution of 2 m surface air temperature (K, contours are −5, −3, −1, 0, 1, 3, 5, 7, 9) averaged over the period 1980–2018 simulated by SOCOLv4 (a, ensemble mean), ERA5.1 reanalysis data (b), and their difference (c, SOCOL4-ERA5.1). Right panels: 2 m surface air temperature trend over the Northern Hemisphere (K/10 years, contours are −0.5, −0.2, 0, 0.2, 0.5, 0.7, 1.0, 1.5) simulated (d) and observed (e) over the Northern Hemisphere. White spots show where the trend is not statistically significant at the 95 % level.

Figure 6Zonal mean zonal wind climatology (m s⁻¹) for DJF (a, c, e) and JJA (b, d, f) simulated with SOCOLv4 (a, b), mean of ERA5.1 and MERRA-2 reanalysis (c, d), and differences between simulated and reanalysis data (e, f). Contour intervals for the difference are −30, −20, −10, −5, −2, 2, 5, 10, 20, 30 m s⁻¹. White areas denote regions where the difference between model and reanalysis data is not statistically significant at the 95 % confidence level calculated with a Student's t test.

Download

Zonal mean climatology of zonal winds and temperatures for boreal winter and summer as simulated by SOCOLv4 are presented in comparison to observations (mean data from reanalysis ERA5.1 and MERRA-2) in Figs. 6 and 7, respectively. Even though the model has an upper limit at 0.01 hPa, we analyze the results only up to a height of 0.1 hPa, since above ERA5.1 data are not reliable and MERRA-2 data are not available. In general, the model reproduces the observations well, but there are deviations which are, however, typical for most coupled models (Bock et al., 2020; Gettelman et al., 2019; Matthes et al., 2020). The tropospheric jets, their position, and strength are reasonably reproduced by the model with biases of ±2–5 m s⁻¹. There is a tendency of the westerlies in the midlatitudes to be insufficiently poleward, which is more pronounced in the Southern Hemisphere and represented as negative and positive deviations around 60 and 40^∘ S, respectively. The differences between the model and observations increase with altitude. In the stratosphere, the simulated summertime easterly jets are of similar strength to those in the observations but have some deviations in their shapes. The deviations between about 90 and 3 hPa in the (sub)tropics are related to the hemispherically symmetric QBO nudging approach, which is applied between 20^∘ N and 20^∘ S (see Sect. 2.4). The wintertime polar night jets show larger deviations in both hemispheres. In the Northern Hemisphere, the modeled polar vortex is significantly weaker than in the observations and is extended equatorward in the lower mesosphere. The Southern Hemisphere vortex is of similar strength to that in reanalysis mean data in the stratosphere and the main model deviations are related to its positioning. The core of the southern polar jet is shifted to 60^∘ and located between 10 and 1 hPa, while in the reanalysis data the vortex core (80 m s⁻¹ isoline) extends higher and shifts towards 40^∘ S. Therefore, the modeled southern polar vortex area is more isolated in the middle stratosphere. Oppositely, in the lower mesosphere, it is extended equatorwards, similar to the northern one. The latter feature might suggest that there is some overestimation in the gravity wave forcing that is dominant in this region (Plumb, 2002; Butchart, 2014).

The zonal mean temperature distribution is generally well reproduced by the model (Fig. 7), but there are also strong biases in the middle atmosphere with respect to the observations that are consistent with the biases in zonal winds. Namely, the wintertime high-latitude middle and upper stratosphere and the low mesosphere are warmer by 2–10 K suggesting weaker vortices and stronger meridional transport and mixing. In the rest of the stratosphere, the difference is negative (about −2 K). This is consistent with a too-fast BDC, which would cause increased adiabatic cooling in the tropics and in the summertime high latitudes. In the lowermost extratropical stratosphere, there are up to −5 K temperature deviations in both hemispheres, which leads to some downward displacement of the extratropical tropopause. The modeled tropospheric temperatures agree very well with observations and only show a warm bias over Antarctica.

It has to be noted that all of the presented model deviations in thermodynamics are long-term issues of the ECHAM-family models and global models in general (Stevens et al., 2013; Mauritsen et al., 2019). In Fig. A1, we present the same differences as in Figs. 6 and 7 but for the pure MPI-ESM CMIP6 run. As can be seen from there, SOCOLv4 and MPI-ESM share the same deviation patterns indicating weaker and warmer polar vortices in the stratosphere. Some differences between SOCOLv4 and MPI-ESM results are mostly related to the stratospheric ozone distribution at midlatitudes and high latitudes, as well as to the distribution of other GHGs except CO₂ that are now three-dimensional, and the QBO nudging in SOCOLv4 in the tropics. Many of these problematic features were shown to be greatly improved using a version of ECHAM with more detailed (95-level) vertical resolution (Schmidt et al., 2013; Stevens et al., 2013; Mauritsen et al., 2019; Matthes et al., 2020), especially the strength of the northern polar vortex and the extratropical lower stratospheric cold biases. The increase of the vertical resolution implies, however, a strong increase in computational time, because of the column physics design. It can be partly compensated by increasing the amount of processing units, but there is also a limit in model's scalability defined by horizontal resolution.

Figure 7Zonal mean temperature climatology for DJF (a, c, e) and JJA (b, d, f) simulated with SOCOLv4 (a, b), mean of ERA5.1 and MERRA-2 reanalysis (c, d), and differences between simulated and reanalysis data (e, f). Contour intervals for the differences are −30, −20, −10, −5, −2, 0, 2, 5, 10, 20, 30, 40 K. White areas denote regions where the difference between model and reanalysis data is not statistically significant at the 95 % confidence level.

Download

3.2.1 Reference data for validation

All datasets that were used for the chemistry validation are listed in Table 1. For the comparison of HCl, H₂O, and N₂O, we used the Global OZone Chemistry And Related trace gas Data records for the Stratosphere database (GOZCARDS, Froidevaux et al., 2015), which is a composite based on the data from eight satellite instruments. The Envisat satellite Michelson Interferometer for Passive Atmospheric Sounding (MIPAS, Funke et al., 2014) instrument data were used for NO_y, while the SCISAT Atmospheric Chemistry Experiment Fourier Transform Spectrometer (ACE-FTS) data climatology is used for NO_y and CH₄. Besides GOZCARDS, the ozone mixing ratio is also compared to the BAyeSian Integrated and Consolidated composite ozone (BASIC, Alsing and Ball, 2019), and the CMIP6 ozone forcing dataset (Checa-Garcia, 2018). The latter is a compilation of CESM1-WACCM and Canadian Middle Atmosphere Model (CMAM) modeling results and is used as a standard forcing for MPI-ESM. Ozone total column is validated against Multi-Sensor Reanalysis version 2 (MSRv2, Van der A et al., 2015a) and Solar Backscatter Ultraviolet Spectral Radiometer version 2 (SBUVv8.6) (McPeters et al., 2013) composites. Stratospheric aerosol sulfur mass is compared to MIPAS (Günther et al., 2018) and High-resolution Infrared Radiation Sounder (HIRS) (Baran and Foot, 1994) observations and a composite recommended by CMIP6 (GloSSACv1.1, Thomason et al., 2018).

It has to be noted that observations themselves are not perfect, especially for the shorter-lived species with pronounced diurnal cycles, low concentrations, and specific regions like high latitudes and the upper troposphere/lower stratosphere (UTLS). Different satellite instruments vary in terms of measurement method, geographical coverage, spatial and temporal sampling and resolution, time period, and retrieval algorithm. All of this has also to be dealt with while constructing the observational composites, which imposes further potential uncertainties related to the homogenization methods applied. The SPARC Data Initiative (SPARC, 2017) conducted an overview of the existing trace gases and aerosol measurements and provided a set of useful recommendations concerning the related uncertainties that we used in the further sections for our validation.

Table 1Data used in this study for chemistry validation.

Download Print Version | Download XLSX

3.2.2 Trace gas climatology

In this section, we discuss the SOCOLv4 performance in reproducing observed distributions of various trace species. Figures 8, 11, and 12 show the data from the model and observations and the relative differences between them. Figures 9 and 13 show mean seasonal cycles of several species that are also compared to observations.

Figure 8Annually averaged zonal mean distributions of CH₄, H₂O, HCl, and N₂O (from top to bottom). First column: SOCOLv4 data. Second column: ACE-FTS and GOZCARDS observational data. Third column: relative difference between the model and observations in %. An observational missing data mask is applied to the modeling results. Note that different species have different averaging periods. White areas are either the missing data or regions where the difference between model and observations is not statistically significant at the 95 % confidence level.

Download

CH₄

Methane (CH₄) is an important source of water vapor (H₂O) in the stratosphere. On average, two H₂O molecules are produced per CH₄ molecule oxidized. Methane mixing ratios drop very fast with height due to increasing availability of hydroxyl radical (OH), the excited atomic oxygen (O(¹D)), and atomic chlorine (Cl). In the lower mesosphere, photolysis by Lyman-α radiation also becomes an important sink for CH₄. Methane is also an important atmospheric tracer for tracking the stratospheric circulation because of its long lifetime. CH₄ mixing ratio isolines in Fig. 8b are lifted upward in the low latitudes by the ascending branch of the BDC and are pushed downward in the middle and high latitudes by the descending branch. Isolines in the midlatitudes are flattened as a result of rapid mixing by breaking planetary waves. The balance between horizontal mixing and the diabatic circulation produces the observed tracer slope and deviations in this distribution can hint at potential problems in transport processes (e.g., Strahan, 2015).

SOCOLv4 nicely reproduces the methane distribution in the lower stratosphere with only about ±5 % deviations from observations. However, above ∼ 10 hPa, the model deviation gradually increases with height to up to 50 % overestimation in the lower mesosphere. The highest overestimation close to the model top in the mesosphere is likely related to the treatment of the Lyman-α line photolysis, while the rest is due to the problems in dynamics. Thus, a too-fast vertical transport in the tropics (Fig. 8a–c) would mean that there is less time for CH₄ loss on the way upwards and therefore more CH₄ arrives in the mesosphere, which can be partly enhanced by underestimation of CH₄ photolysis. Dynamical lifetime got smaller with respect to chemical lifetime and the tracer became more vertically mixed. The increased concentration in the upper stratosphere is then also distributed polewards and downwards in the middle and high latitudes. Difference in skewness of the isolines suggests an overestimation of midlatitude mixing in the upper stratosphere/lower mesosphere, which is consistent with analysis of winds showing that the vortices are weaker at those altitudes (Fig. 6). Middle and lower stratospheric midlatitudes, however, show some hemispheric asymmetry with a negative bias of 5 % in the SH and a positive bias of 5 %–10 % in the NH. Figure 9a–b show the CH₄ seasonal cycle at 20 hPa and suggests that this anomaly is likely related to the deficiencies in horizontal mixing. Namely, while the model reproduces very well the tropical values (1.4 ppmv isoline), the midlatitude values (1 ppmv isoline) are excessively shifted poleward in the NH and slightly equatorward in the SH. This is also generally consistent with Fig. 6, which shows at these altitudes a weaker NH polar vortex and the SH polar vortex that is more confined to the pole.

Figure 9CH₄, HCl, and NO_y zonal mean seasonal cycles at 20, 70, and 6 hPa, respectively. First column: SOCOLv4 data. Second column: ACE-FTS, GOZCARDS, and MIPAS observational data. Two consecutive cycles are shown.

Download

H₂O

Water vapor (H₂O) provides a source for hydroxyl radicals (HO_x) in the stratosphere. Hydroxyl catalytic cycles dominate the ozone destruction above the stratopause and are also very important in the lower and middle stratosphere. The two main sources of H₂O in the stratosphere are methane oxidation and the transport from the H₂O-rich troposphere. Transport from the troposphere is important for the lower stratosphere and is highly dependent on the cold point temperature at the tropical tropopause and vertical transport in the tropics. As seen from Fig. 8d–f, SOCOLv4 agrees very well with observations around the stratopause. In the lower stratosphere, the model is 5 %–10 % too moist, which is due to some overestimation of H₂O at the entry point. In the mesosphere, the model underestimates H₂O by up to 30 % (or 1 ppmv) around the model top. This partly results from the underestimated destruction of methane (about 0.2 ppmv) but mostly comes from the overestimated water vapor photolysis, which was also previously reported for SOCOL in the stand-alone photolysis codes intercomparison study (Sukhodolov et al., 2016; Karagodin-Doyennel et al., 2021). The water vapor sink (dehydration) by PSCs in the south polar lowermost stratosphere is well captured by the model.

Water vapor is a relatively long-lived species with a highly pronounced seasonality in the lower tropical stratosphere. This seasonality is related to tropospheric wave forcing, which is strongest in northern winter and weakest in northern summer, leading to enhanced upwelling and colder temperatures at the tropical tropopause during boreal winter, and vice versa. In combination with the tropical pipe and the subtropical transport barrier, this makes H₂O useful as a tracer of transport processes in the tropical stratosphere (Mote et al., 1998). Figure 10 shows the water vapor “tape recorder” signals obtained from GOZCARDS and SOCOLv4 to illustrate the model performance in terms of vertical transport in the tropical low and middle stratosphere. SOCOLv4 reproduces the seasonal cycle of water vapor at the tropopause level well. However, it also shows too-fast vertical propagation of the water vapor tape recorder in the observations similar to its predecessor versions (Stenke et al., 2013), which is a common peculiarity for chemistry–climate models (e.g., Dietmüller et al., 2018). This issue is related to the problems in horizontal and vertical mixing and diffusion and also to representation of the residual circulation in general, all of which are largely dependent on resolution. The overestimation of H₂O of about 0.5 ppmv throughout the middle and lower stratosphere (see also Fig. 7) is mostly due to the higher entry value during summer but can also be partly related to the overestimation of mixing with higher altitudes.

Figure 10Annual cycle of water vapor mixing ratio (ppmv) averaged between 15^∘ N and 15^∘ S and over 1991–2012 for the GOZCARDS satellite composite (a) and SOCOLv4 (b). Two consecutive cycles are shown. The white line in both panels indicates the approximate slope of the GOZCARDS tape recorder signal.

Download

HCl

HCl is a chlorine reservoir that is widely used for tracking the overall chlorine abundance in the stratosphere, which also plays an important role in polar ozone chemistry. SOCOLv4 reproduces the HCl climatology pretty well with deviations mostly on the order of ±5 %–10 % throughout the stratosphere (Fig. 8g–i). The most pronounced deviation is the ∼ 5 % underestimation by the model at the stratopause and the tropical mid-stratosphere. Figure 9c–d illustrate the seasonal cycle at 70 hPa, demonstrating a rather good representation of the southern high-latitude variations with some underestimation of the PSC activation effects, the reasons for which are yet unclear. A limiting factor here can be the amount of available ClONO₂ (Nakajima et al., 2016), whose observations are, however, rather sparse for these conditions. Note that HCl observational uncertainties are also the largest in the SH polar vortex (SPARC, 2017). In the northern polar regions, SOCOLv4 reproduces the semi-annual variability seen in observations but overestimates the wintertime decrease and underestimates the summertime one, which might be related to the weaker northern vortex that leads to higher temperature and lower frequency of PSC formation.

NO_y

Figure 11 shows the odd nitrogen (NO_y = NO + NO₂ + 2N₂O₅ + HNO₃) climatology from the model and its difference to the ACE-FTS and MIPAS instruments data. Odd nitrogen provides another important sink of ozone in the stratosphere and also negatively contributes to and is affected by other catalytic cycles of ozone destruction by deactivating odd hydrogen, chlorine, and bromine in the reservoir species HNO₃, ClONO₂, and BrONO₂, respectively. The main source of NO_y in the middle atmosphere is oxidation and photolysis of nitrous oxide (N₂O). The lowermost stratosphere is additionally highly affected by the mixing with the troposphere and therefore the NO_y budget also depends on the free tropospheric sources of odd nitrogen such as aircraft engines and lightning (Grewe, 2009). In the mesosphere, where odd nitrogen is mostly present in the form of nitric oxide (NO), the most important sources become the in situ N₂ ionization by solar protons and middle energy electrons (MEEs) and the transport from the thermosphere, where it is produced by auroral electrons and N₂ photolysis. Compared to ACE-FTS and MIPAS, which generally agree with each other, SOCOLv4 overestimates by ∼ 10 % the maximum concentrations in the middle–upper stratosphere, the region where the odd nitrogen catalytic cycle of ozone destruction is the most effective. In the upper stratosphere and the mesosphere, the model also overestimates NO_y everywhere except polar mesospheric regions. The overestimation increases gradually with height, reaching 80 % in the tropical mesosphere, where the absolute values, however, are rather small. As there is almost no N₂O left at these altitudes, the overestimation likely comes from the underestimated sink through the NO photolysis. Despite the strong percentage deviation in the tropical mesosphere, the model underestimates NO_y in its polar regions, though the differences are almost insignificant due to a large interannual variability. The polar bias is most likely related to an underestimated contribution of either MEEs or the thermospheric NO downward flux. Mironova et al. (2019) compared CMIP6 MEE ionization effects to calculations based on polar balloon measurements and found that the related odd nitrogen signal can be up to 100 % too low when CMIP6 forcing was used. High-latitude and midlatitude upper stratospheric NO_y concentrations are then affected by both the positive bias in the low latitudes and the negative bias in the polar high latitudes.

Figure 11NO_y zonal mean annual mean climatology (ppbv) calculated by (a) the SOCOLv4 model and (b) provided by MIPAS. (c, d) Relative difference between the model and the ACE-FTS and MIPAS observations. White areas are either the missing data or regions where the difference between model and observations is not statistically significant at the 95 % confidence level.

Download

NO_y observational uncertainty is the largest in the lower stratosphere (±30 %, SPARC, 2017), however both instruments suggest that the model has a pronounced low bias there. This can be explained by the insufficient production from N₂O oxidation, which shows an underestimation in the middle and upper stratosphere (Fig. 7). This underestimation is very similar to that in methane, which suggests that both are probably related to the overestimation of the vertical transport. Another potential reason could be overestimated sinks to reservoir species; however, modeled HNO₃ and ClONO₂ (Fig. A2) agree reasonably well (±10 %) with observations in the areas of their maximum concentrations in the lower stratospheric midlatitudes. The middle stratospheric overestimation of HNO₃ is likely related to the overestimated abundance of H₂O and hence HO_x at these altitudes and the underestimated sink through photolysis (Sukhodolov et al., 2016). The tropical UTLS region shows large underestimation of both NO_y and HNO₃, which suggests that either the free-tropospheric sink of odd nitrogen through HNO₃ washout or the source through production by lightning activity is biased. Overall, NO_y representation in the model is defined by an underestimated sink in the upper stratosphere and the mesosphere, and the underestimated production in the lower stratosphere, which results in the slight vertical shift of the maximum NO_y concentrations upwards.

Modeled NO_y seasonal cycle at 6 hPa is presented in Fig. 9e–f. The model shows generally higher values than MIPAS due to NO_y burden being biased high, as discussed before. However, the seasonality of the variations is reasonably well reproduced at this altitude. The high-latitude and midlatitude seasonality is controlled by the vertical motions, which either bring the NO_y-low air from above during the wintertime or the NO_y-high air from below in the summertime. This transport is modulated by horizontal mixing, which can cause some time shifts in the appearance of minima, like seen in the NH winter modeled results, where the contribution of the vertical transport is largely diminished by the mixing with midlatitudes. In the southern vortex, this is better reproduced by SOCOLv4 for the upper stratosphere; however, in the middle and lower stratosphere (not shown), there are traces of the modeled southern vortex being too isolated, which is consistent with the analysis of CH₄ and N₂O. Boreal winter minimum in the tropics appears due to increased upwelling, and thus transport from below, which is well caught by the model.

Figure 12Ozone zonal mean annual mean climatology (1985–2012) calculated by the SOCOLv4 model (ppmv) (a) and provided by the GOZCARDS data (b). (c, d) Relative difference between the model and the observational composites BASIC and GOZCARDS. White areas are either the missing data or regions where the difference between model and observations is not statistically significant at the 95 % confidence level.

Download

O₃

Figure 12 represents the ozone mixing ratio climatology calculated by the model versus two different ozone datasets. Both observational datasets (BASIC and GOZCARDS) agree very well with each other. The difference plots demonstrate that SOCOLv4 has a very good representation of ozone in the tropics and subtropics throughout the stratosphere, except some overestimation in its lowermost part. The lower stratospheric overestimation spreads further to midlatitudes and reaches values of more than 20 %. This overestimation is a common feature in many CMIP6 models (Keeble et al., 2021 – see their Fig. 4). The lower stratosphere is a complex region with competing transport, mixing, and chemistry processes, variations of the tropopause, and exchange with the troposphere (Holton et al., 1995), which makes it difficult to identify exact factors responsible for these biases. From the chemistry side, there is no clear evidence from our analysis that can explain this overestimation. The low NO_y model bias (Fig. 11) would also contribute negatively to the ozone burden at these altitudes through reduced production by oxidation of carbon monoxide, methane, and higher hydrocarbons, similarly to the troposphere (Gauss et al., 2006; Grewe, 2009). Therefore, it is most likely that the high ozone bias is related to the photolysis problems in SOCOLv4. Sukhodolov et al. (2016) analyzed the stand-alone performance of the lookup-table photolysis scheme applied in SOCOL and found that the ozone photolysis in both channels is underestimated especially in the lower stratosphere. Although that study also showed the underestimated molecular oxygen photolysis at the same altitudes, which could partly compensate for the ozone photolysis underestimation, it appears that the latter still dominates. Comparison to the results of other codes in Sukhodolov et al. (2016) showed that the lower stratospheric problems in photolysis are probably related to the neglected Rayleigh scattering effect and temperature dependence of absorption cross sections and quantum yields. The applied photolysis corrections, as discussed in the MEZON module description, seem to improve ozone throughout the stratosphere; however, a more qualitative treatment of photolysis is already planned for the future. Insufficient ozone photolysis also might affect HO_x formation via O(¹D) + H₂O and related catalytic ozone destruction, thus intensifying the ozone high bias. Other factors that could contribute to the lower stratospheric overestimation in ozone are the ∼ 2 K underestimation in local temperatures (Fig. 6), which would slow down the ozone destroying chemistry, and overestimated tropospheric ozone. In addition, SOCOLv4 so far does not include iodine chemistry, while some recent studies suggested that it can be responsible for ∼ 30 % of the halogen-induced ozone loss in the lowermost stratosphere (e.g., Koenig et al., 2020).

In the upper stratosphere at midlatitudes, SOCOLv4 underestimates ozone by up to −20 %. This feature has been present in SOCOL models since the first version (Egorova et al., 2005) and is related to the transport issues, caused by the insufficient vertical resolution affecting the wave propagation and breaking in the stratosphere/mesosphere. Though SOCOLv4 now has 47 vertical levels compared to 39 in previous versions, all additional levels are located in the troposphere, while the number of levels in the middle atmosphere remains almost the same. As was discussed earlier, SOCOLv4 and its atmospheric base (ECHAM6) show stronger residual circulation that causes stronger vertical transport in the tropics and also stronger meridional and downward transport in the winter hemisphere. The increased BDC, however, also implies a stronger wave forcing for the polar vortex making it weaker and blurring the midlatitude transport barriers (Butchart, 2014). An associated increase of stirring, which dominates over the horizontal transport by the meridional circulation in these conditions (Plumb, 2007), would then lead to a decrease in the meridional gradient of tracers. As depicted in Fig. 6, the modeled polar vortex strength in both hemispheres is notably underestimated in the uppermost stratosphere and the lower mesosphere.

Figure 13 shows that the midlatitude ozone underestimation at 1 hPa appears in the winter months and that it happens together with an underestimation in the polar-most latitudes due to the more intense downward transport from the mesosphere. These results are consistent with those seen from the analysis of methane. Part of the model biases here could be related to the 30 % NO_y overestimation (Fig. 11) that would intensify the ozone destruction. In the lower levels (10 and 70 hPa), the model shows a generally good agreement with GOZCARDS, with some deviations in timing and magnitude of extremes that can also be partly attributed to the deficiencies in transport. Thus, a too-deep ozone drop at 10 hPa in the southern polar vortex is consistent with the too-low methane values under these conditions and can be explained by the insufficient mixing. Seasonality of the tropical values at all latitudes agrees very well with observations, showing a decrease of ozone in the lower stratosphere during boreal winter and a pronounced semi-annual oscillation in the middle and upper stratosphere.

Although the observations are very uncertain in the polar night regions, both the ozone mixing ratio (Fig. 13) and the total ozone measurements (Fig. 14) suggest that in SOCOLv4 the size and intensity of the Antarctic ozone hole is overestimated. This can be related to either a too-strong vortex isolation, as discussed before, or to problems in the polar heterogeneous chemistry. Since the low bias is also apparent in the higher levels, another potential reason could be the poleward dislocation of the southern polar night jet in the upper and middle stratosphere (Fig. 6), which would imply less upper stratospheric ozone production inside the vortex (McConnel and Jin, 2008) and thus less ozone available for the downward transport. The downward transport itself might also be biased high, as illustrated by small negative anomalies of N₂O and CH₄ in Fig. 8, which would negatively contribute to ozone by bringing too much chlorine from above.

A comparison of SOCOLv4 ozone to the prescribed CMIP6 ozone used by MPI-ESM (Fig. A3) is useful for understanding the differences in temperature and winds between SOCOLv4 and MPI-ESM in Fig. A1. The comparison shows that SOCOLv4 has slightly more ozone in the tropical and subtropical upper stratosphere but less ozone in the wintertime high latitudes. SOCOLv4 also has higher ozone in the mesosphere (up to 70 %) and in the troposphere (up to 20 %). Upper and mid-stratospheric differences in ozone contribute to a slightly stronger meridional temperature gradient due to less heating in the high latitudes and more heating in the tropics. This makes the polar vortices, which are too weak in MPI-ESM, slightly stronger in SOCOLv4. It is the most pronounced for the SH polar vortex, which is weaker than in observations throughout the whole stratosphere in MPI-ESM (Fig. A3d), while in SOCOLv4 (Fig. 6f) the negative anomaly vanished in the middle stratosphere.

Figure 13The 1, 10, and 70 hPa ozone zonal mean seasonal cycles (ppmv). First column: SOCOLv4 data. Second column: GOZCARDS observational data. Two consecutive cycles are shown.

Download

3.2.3 Total ozone climatology

Total ozone column (TOC) represents an aggregate of all features and model biases discussed above, as well as their compensating effects. Figures 14 and 15 show that SOCOLv4 has a very good representation of tropical and subtropical TOC compared to MSRv2 and SBUV data, which is consistent with the results of Figs. 12 and 13. The model accurately simulates both the seasonal cycle (Fig. 14) and the annual mean climatologies (Fig. 15). Mostly due to the corrected photolysis rates, SOCOLv4's TOC is about 20 DU lower in the tropics and subtropics than in the previous model version, which results in better agreement with observations. Compared to MSRv2 data, SOCOLv4 shows almost perfect annual mean TOC in the tropics with some regional biases on the order of ±10 DU (see also Fig. A4 for TOC anomalies). Midlatitudes are also generally well captured, with the most pronounced overestimation of up to 30 DU being the spring seasons in the Northern Hemisphere. Note also the ±10 DU disagreement between SBUV and MSRv2 in the midlatitudes.

Figure 14Seasonal cycle of total ozone simulated by SOCOLv4 (a) against SOCOLv3 (b), MSRv2 (c), and SBUV observations (d) averaged for the 1980–2014 period.

Download

Major model bias appears in the polar latitudes under wintertime conditions. It has to be mentioned that the polar night observations are subject to substantial uncertainties. SBUV data give no information for the polar night conditions at all, while MSRv2 provides almost full spatial and temporal coverage by assimilating a large variety of data including the ground-based observations into a chemistry-transport model (Van der A et al., 2015a). This model, however, exploits a simplified scheme of stratospheric chemistry by Cariolle and Teyssèdre (2007) that tends to overestimate TOC under the ozone hole conditions by about 20 DU. However, taking this uncertainty into account while also using the direct satellite information by SBUV and GOZCARDS (Fig. 12) at the edge of the polar vortices, SOCOLv4 still shows some pronounced biases. Namely, the polar ozone depletion in the SH starts right after the vortex formation (Fig. 14a), while the observations and the current state of knowledge suggests that it starts later under the springtime conditions (Fig. 14c). This feature is fully inherited from the previous model version (Fig. 14b) and needs further detailed investigation. Likely reasons, as discussed earlier, could be related to the transport issues, the details of polar heterogeneous chemistry, or also the treatment of photolysis under high solar zenith angles. In the Northern Hemisphere, Figs. 6 and 7 suggest a clear underestimation of the NH vortex strength in the model, which can explain the difference with respect to observations. The observations suggest some local ozone decrease in the middle of the winter, which is then followed by a fast increase in the springtime caused by inflow from the midlatitudes. The vortex in the model is weaker, meaning that polar air masses are less isolated and can get enriched with ozone from midlatitudes throughout the wintertime, and therefore the modeled field shows a more gradual increase instead of an abrupt rise. This stronger mixing of polar vortex air is also responsible for the NH results in Fig. 15, where the model correctly reproduces the shape of the high ozone values but also shows some overestimation in the modeled annual mean total ozone.

Figure 15Global distribution of annual mean total ozone for 1980–2018 as (a) simulated by SOCOLv4 and (b) from the MSRv2 data record in DU. (c) Difference of panels (a)–(b) in DU with contour lines −40, −30, −20, −10, −5, 0, 5, 10, 20, 30 DU.

3.2.4 Ozone evolution

The upper panel of Fig. 16 shows the evolution of the total ozone column (TOC). Near-global TOC was calculated for this comparison, because SBUV has no data over polar night regions, and the reliability of MSRv2 there is also low. SOCOLv3 and SOCOLv4 are presented as three separate ensemble members to illustrate the uncertainty due to the internal variability. SOCOLv3 data are taken from the RefC1 runs of CCMI-1, which assume prescribed historical sea surface temperatures and sea ice coverage. Using prescribed sea surface temperatures (SSTs) in v3 explains why the variability among ensemble members is smaller in v3 than in v4, where the ocean is now interactive, and the ENSO variability is stochastic rather than fixed (Fig. 4). It is clearly seen that SOCOLv3 has a large overestimation in TOC compared to other time series. As was discussed above, this is mostly related to the problems in the photolysis scheme that are now manually corrected. SOCOLv4 and MSRv2 agree well in terms of the mean values, while the mean TOC from SBUV is lower by about 5–7 DU with respect to MSRv2. Generally, the slight overestimation of global TOC is also seen in many CMIP6 models, and SOCOL is well within the spread across CMIP6 models (Keeble et al., 2021, see their Fig. 6a). The CMIP6 ozone composite, instead, shows slightly lower TOC than SOCOLv4 and MSRv2 in Fig. 16a.

Figure 16(a) Near-global (55^∘ S–55^∘ N) mean ozone evolution simulated with SOCOLv4 (black), MSR2 (red), SBUV (dark blue), SOCOLv3 (light blue), and CMIP6 (green). (b) The same but anomalies with respect to the mean over the whole period.

Download

For easier comparison of the TOC variability, we also show anomalies in Fig. 16b. The four decades shown are characterized by contributions from anthropogenic and natural forcings with varying importance during different periods. From 1980 to the mid-1990s, there is a strong hODS-induced ozone decline. This negative tendency is further enhanced by the two strong volcanic eruptions of El Chichón in 1982 and Pinatubo in 1991, which are then followed by several years of recovery. Since the mid-1990s, following a recovery after Pinatubo, ozone starts to increase further, induced by the leveling-out of the hODS emissions as a consequence of the Montreal Protocol (Egorova et al., 2013; Chipperfield et al., 2015). Since the 2000s, there have been no major eruptions, and the ozone evolution is mostly determined by the decreasing halogen load and the global warming effects, though there are still many related uncertainties and other important factors in specific regions, especially in the lower stratosphere (Ball et al., 2018; Stone et al., 2018). Upper stratospheric ozone was shown to have the most pronounced recovery (Petropavlovskikh et al., 2019; Chipperfield et al., 2018) due to being the most sensitive to the chlorine catalytic cycle of ozone destruction. Furthermore, decreasing stratospheric temperatures also led to some deceleration of ozone destruction cycles. The contribution of these positive changes to the total column is partly supplemented by an increase of the tropospheric column ozone at a rate of about 6 %–7 % per decade resulting from the continuous increase of the surface anthropogenic emissions of ozone precursors like CH₄, CO, NO_x, and VOC (Ball et al., 2018). The recovery in the middle and lower stratosphere in the tropical and extratropical areas is still unclear and observations even suggest negative trends in the lowermost stratosphere, which are, however, barely significant due to a large dynamical variability of this region (Ball et al. 2018; 2019; Petropavlovskikh et al., 2019). There are indications that these negative trends are dynamically driven (Orbe et al, 2020) and can be related to acceleration of the BDC due to rising GHGs; however, models mostly fail to reproduce the observations in detail (Ball et al., 2020).

In Fig. 16b, we show that SOCOLv3 had a reasonable TOC evolution under volcanically quiescent conditions, while the volcanic effects were strongly overestimated. Unlike SOCOLv4, which uses an interactive sulfate aerosol scheme, SOCOLv3 used prescribed fields recommended for CCMI that seem to be biased, probably because of the uncertain observations and also the assumptions made during the calculation of aerosol mass from extinction observations. The interactive aerosol scheme clearly improves the performance of SOCOLv4 in terms of volcanic effects on ozone. Sukhodolov et al. (2018) showed though that the observed hemispheric asymmetry of volcanic signals in ozone are still difficult to properly reproduce with a free-running model. Overall, the general behavior of the TOC evolution is very well captured by the new model both in terms of the mean values of the near-global field and the variability.

In Fig. 16, we looked into some drivers of the total column evolution. Depletion and further recovery of ozone in the middle stratosphere are well captured by the model (Fig. 17a) due to the good representation of changes in stratospheric chlorine (Fig. 17e) and the continuous cooling of the stratosphere (Fig. 17b). The lower stratospheric (60 hPa) temperature and ozone observations show a much larger uncertainty but still somewhat pronounced negative trends that are also well captured by the model. A thorough analysis of lower stratospheric processes requires comprehensive statistical tools, and a set of sensitivity runs, which is already planned for future model applications.

Figure 17(a, c) Near-global (55^∘ S–55^∘ N) annual mean ozone evolution at 60 and 10 hPa compared to the BASIC composite. (b, d) Tropical annual mean evolution of temperature at 60 and 10 hPa compared to ERA5. (e) Tropical monthly mean (smoothed with the 3-month running mean) evolution of HCl at 10 hPa compared to the GOZCARDS composite. Black lines: three model ensemble members. Orange lines: observations.

Download

Figure 18Global total stratospheric sulfur mass in sulfate aerosols (a) and its zonal mean values (b–d) simulated with SOCOLv4 against HIRS, CMIP6, and MIPAS. Data are plotted in a logarithmic scale so that both large and small eruptions can be distinguished.

Download

3.2.5 Stratospheric sulfur

Figure 18 shows the total stratospheric sulfur load as well as the zonal mean stratospheric column of liquid H₂SO₄ for SOCOLv4 against CMIP6, MIPAS, and HIRS. The tropopause that was applied to calculate the stratospheric burden for the CMIP6 and MIPAS datasets was calculated from the ERA-Interim temperature data using the WMO definition of the tropopause as the lowest level at which the temperature decreases less than 2 K km⁻¹, while for SOCOLv4 the modeled tropopause is used. The peaks in aerosol originate from volcanic eruptions with the largest being the El Chichón (1982) and Pinatubo (1991) eruptions in the earlier time period and Kasatochi (2008), Sarychev (2009), Nabro (2011), and Calbuco (2015), as well as some others in the more recent past.

Generally, the aerosol evolution in the stratosphere is depicted quite well by the model in terms of the total aerosol amount, lifetime, and meridional transport. There are, however, a few differences that need to be addressed. Prior to and after the 1982 El Chichón eruption, there are some significant discrepancies between SOCOLv4 and CMIP6. In addition, the peak of stratospheric burden caused by this eruption (see Fig. 18a) is slightly shifted in time in comparison to the SOCOLv4 simulations. In Fig. 18c, the delayed peak of CMIP6 is seen mostly in the Northern Hemisphere and the tropical aerosol levels, where the eruption took place, are lower than in SOCOLv4. The El Chichón eruption occurred between the activities of the Stratospheric Aerosol and Gas Experiment (SAGE) and SAGE II. Data from this time are a composite of Stratospheric Aerosol Measurement II (SAM II) data in the extratropics, as well as several airborne lidar missions at low latitudes, which delivered data often many months apart, leading to very sparse data coverage. The gaps were then filled by means of interpolation (Thomason et al., 2018). It is, therefore, difficult to make a conclusive statement on the SOCOLv4 performance of this time due to large uncertainties in observations. The next larger event, Nevado del Ruiz in 1985, already occurred during the SAGE II era, which is deemed more reliable by Thomason et al. (2018).

The 1991 Mt. Pinatubo eruption has been discussed several times, since observations show a plateaued value for several months after the eruption instead of a more distinct peak as models like SOCOLv4 suggest (e.g., Dhomse et al., 2020). Something similar has not been observed since then. However, no other well-observed eruption was the size of the 1991 Mt. Pinatubo one. A source for measurement artifacts during this time is the increase in opacity of the atmosphere due to large amounts of aerosol resulting from the eruption. From the modeling side, a potential deficiency could be insufficient vertical resolution of the stratosphere leading to increased vertical diffusion and faster meridional transport of aerosols. Niemeier and Schmidt (2017) used a model with a similar dynamical core, ECHAM-HAM, with 39 and 90 vertical levels. They showed that the version with a higher vertical resolution maintains aerosols longer in the tropical reservoir, which increases their lifetime after low-latitude emission events. The post-plateau e-folding time looks very similar for CMIP6 and SOCOLv4, but in the CMIP6 data there is another small increase after the main Mt. Pinatubo plateau without there having been any other major eruptions. Again, this could be due to a change in the instruments used for the construction of this dataset.

The volcanically quiescent period in the late 1990s and early 2000s and the minor volcanic activity later, which is also better covered by observations, are well reproduced by the model. MIPAS observations are added here, as they are not included in the CMIP6 composite and thus represent an independent source of information. MIPAS and CMIP6 agree with each other quite well in terms of elevated values after eruptions, while the background level is slightly lower in MIPAS. Background levels in SOCOLv4 are located between CMIP6 and MIPAS and there is some spread over the ensemble members, which is caused by the variable tropopause in the model combined with vertically prescribed SO₂ emission profiles. The model also captures the zonal mean distribution of elevated aerosol values, reproducing some equatorward transport after high-latitude eruptions of Kasatochi, Sarychev, and Calbuco, as well as dominant northward transport after the equatorial eruption of Nabro. It should be noted that our modeling results here are largely defined by the emission database that we used (Carn et al., 2016). Other volcanic emission databases report quite different estimates for each volcanic event, both in terms of emitted amount of SO₂ and its vertical distribution (Timmreck et al., 2018).

Overall, we can conclude that SOCOLv4 nicely reproduces background aerosol levels and minor volcanic activity compared to observations, while potentially underestimating the aerosol lifetime after major eruptions of El Chichón, Nevado del Ruiz, and Pinatubo. The latter can be attributed to a number of factors. As mentioned above, the data from earlier periods are less reliable and had a lower temporal coverage (Thomason et al., 2018), which influences the forcing data used in the model, as well as the observational composites for its validation. Furthermore, the physical processes influencing aerosol concentrations after these larger eruptions differ from the smaller but better observed ones that happened more recently. Since more SO₂ is injected into the stratosphere at once in a larger event, OH radicals are depleted more quickly, which can cause a delay in sulfate aerosol formation. Nevertheless, all SO₂ is eventually (in 2–3 months' time) converted to sulfate aerosol and in higher concentrations the particles will collide and grow in size. These larger aerosols are subjected more to the force of gravity and more easily sediment to the troposphere (Timmreck et al., 2009), while the long-term decay of aerosols after smaller eruptions is more controlled by the BDC transport (Günther et al., 2018). Given that all major eruptions were tropical, a bias in vertical numerical diffusion due to insufficient vertical resolution could further intensify the sedimentation to the lower stratosphere and the troposphere and also reduce the confinement of aerosols in the tropical reservoir via increased leakage through the shallow branch of the BDC.