The chemistry–climate model ECHAM6.3-HAM2.3-MOZ1.0

The chemistry–climate model ECHAM-HAMMOZ contains a detailed representation of tropospheric and stratospheric reactive chemistry and state-of-the-art parameterizations of aerosols using either a modal scheme (M7) or a bin scheme (SALSA). This article describes and evaluates the model version ECHAM6.3-HAM2.3-MOZ1.0 with a focus on the tropospheric gas-phase chemistry. A 10-year model simulation was performed to test the stability of the model and provide data for its evaluation. The comparison to observations concentrates on the year 2008 and includes total column observations of ozone and CO from IASI and OMI, Aura MLS observations of temperature, HNO₃, ClO, and O₃ for the evaluation of polar stratospheric processes, an ozonesonde climatology, surface ozone observations from the TOAR database, and surface CO data from the Global Atmosphere Watch network. Global budgets of ozone, OH, NO_x, aerosols, clouds, and radiation are analyzed and compared to the literature. ECHAM-HAMMOZ performs well in many aspects. However, in the base simulation, lightning NO_x emissions are very low, and the impact of the heterogeneous reaction of HNO₃ on dust and sea salt aerosol is too strong. Sensitivity simulations with increased lightning NO_x or modified heterogeneous chemistry deteriorate the comparison with observations and yield excessively large ozone budget terms and too much OH. We hypothesize that this is an impact of potential issues with tropical convection in the ECHAM model.

1 Introduction

Global chemistry–climate models have become indispensable tools for the investigation of interactions between atmospheric chemistry and various aspects of the physical and biogeochemical climate system. In recent years, several coupled models have been developed with varying levels of interaction between Earth system compartments and varying details in their representation of chemical and physical processes (Young et al., 2013; Morgenstern et al., 2017; Young et al., 2018).

Here, we describe and evaluate a new chemistry–climate model based on the general circulation model ECHAM6.3 (Stevens et al., 2013), the Hamburg Aerosol Model (HAM) version 2.3 (Tegen et al., 2018; Stier et al., 2005; Zhang et al., 2012), and the gas-phase tropospheric and stratospheric module MOZ1.0.

ECHAM6.3-HAM2.3-MOZ1.0 (henceforth ECHAM-HAMMOZ) can be run in different configurations: (1) using prescribed fields of surface pressure, divergence, vorticity, and temperature and applying a relaxation technique with time-varying weights (“nudging”); (2) constraining only sea surface temperatures and sea ice concentrations (“AMIP mode”); or (3) fully coupled with ocean and sea ice models. In this study we concentrate on simulations of type 1 as these allow for a more detailed evaluation of the model with observational data and because most current applications of ECHAM-HAMMOZ use this mode. For a discussion on the potential differences in the different configurations the reader is referred to Lamarque et al. (2012).

Earlier versions of ECHAM-HAMMOZ have been used successfully to analyze the impact of heterogeneous reactions on tropospheric ozone chemistry (Pozzoli et al., 2008a) and on aerosol composition (Pozzoli et al., 2008b) over the North Pacific, the influence of African emissions on regional and global tropospheric ozone (Aghedo et al., 2007), the impact of continental pollution outflow on the chemical tendencies of ozone (Auvray et al., 2007), and the impact of Asian aerosol and trace gas emissions on the Asian monsoon (Fadnavis et al., 2013, 2014, 2015). A 25-year reanalysis with ECHAM-HAMMOZ was performed by Pozzoli et al. (2011). In addition, several studies were performed with the aerosol climate model ECHAM-HAM, which uses trace gas climatologies from ECHAM-HAMMOZ to constrain aerosol nucleation (e.g., Jiao et al., 2014; Neubauer et al., 2014; Stanelle et al., 2014; Ghan et al., 2016; Zhang et al., 2016). The tropospheric chemistry–climate model ECHAM5-MOZ also participated in the first multi-model intercomparison study of the Task Force Hemispheric Transport of Air Pollutants (TFHTAP) (Dentener et al., 2006a; Stevenson et al., 2006).

This article intends to provide a thorough description of the chemistry component of the ECHAM-HAMMOZ model with special focus on tropospheric reactive gases. Stratospheric chemistry is briefly discussed as well, while for more detailed discussions of the performance of the physical climate model ECHAM6.3 the reader is referred to Stevens et al. (2013). More information on the aerosol schemes HAM-M7 and HAM-SALSA and their evaluation can be found in Stier et al. (2005), Zhang et al. (2012), Neubauer et al. (2014), Tegen et al. (2018), and Kokkola et al. (2018).

This article first provides general descriptions of the ECHAM6.3, HAM2.3, and MOZ1.0 components (Sect. 2) before the gas-phase chemistry parameterizations are discussed in more detail (Sect. 3). Section 4 provides an overview of the simulations performed for this paper. Section 5 presents simulation results and comparisons with observations and other independent model simulations. In Sect. 6 we analyze the global budgets of ozone, OH, NO_x, aerosols, clouds, and radiation. Section 7 contains conclusions.

2 Model description

2.1 ECHAM6.3

ECHAM6, subversion 3, is the sixth-generation general circulation model from the Max Planck Institute for Meteorology in Hamburg, Germany (Stevens et al., 2013). The model uses a spectral dynamical core to calculate temperature, surface pressure, vorticity, and divergence. Diabatic processes such as convection, diffusion, turbulence, and gravity waves are calculated on an associated Gaussian grid. The vertical discretization is a hybrid sigma–pressure coordinate system.

Transport of scalar quantities is performed with the flux-form semi-Lagrangian scheme of Lin and Rood (1996). Turbulent mixing adopts an eddy diffusivity and viscosity approach following Brinkop and Roeckner (1995), and 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). Cloud water and cloud ice are treated prognostically according to Lohmann and Roeckner (1996). In the base model version, the cloud droplet number concentration is parameterized as a function of altitude with higher values over land than over the ocean. In contrast, ECHAM-HAMMOZ explicitly calculates cloud droplet number concentration as a function of aerosol activation (see below). Gravity waves are generated from a subgrid orography scheme (Lott, 1999) and as Doppler waves following Hines (1997a, b), and they are treated according to the formulation of Palmer et al. (1986) and Miller et al. (1989). Radiative transfer calculations are done with the two-stream method of RRTM-G (Iacono et al., 2008). 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 ECHAM-HAMMOZ uses the prognostic tracer concentrations of aerosol and the following gases to specify absorption and scattering: CO₂, CH₄, N₂O, CFC11, CFC12, O₂, and O₃. 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).

Land surface processes are modeled with JSBACH (Reick et al., 2013), which uses a tiling approach with 12 plant functional types and two types of bare surface. The soil hydrology and temperatures are modeled by a five-layer scheme (Hagemann and Stacke, 2015), which constitutes an update from the description provided by Stevens et al. (2013).

2.2 HAM2.3

The Hamburg Aerosol Model (HAM) consists of parameterizations of all relevant aerosol processes including emissions, nucleation, condensation, coagulation, cloud activation, dry deposition, wet deposition, and sedimentation. HAM solves prognostic equations for sulfate, black carbon, particulate organic matter, sea salt, and mineral dust aerosol. Two different representations of aerosol microphysics are available based on the modal scheme M7 (Vignati et al., 2004; Stier et al., 2005) or on the Sectional Aerosol module for Large Scale Applications (“SALSA”: Kokkola et al., 2008; Bergman et al., 2012). These microphysical packages solve for the tendencies of nucleation, condensation, coagulation, and hydration. Since previous versions of ECHAM-HAM and ECHAM-HAMMOZ only included the modal approach of M7, the aerosol processes emissions, wet and dry removal, particle-phase chemistry, and radiative properties were generalized to also function with the sectional approach. The simulations described in this paper were performed with the M7 scheme. M7 represents aerosol sizes as seven modes, i.e., four soluble and three insoluble modes with fixed standard deviation, but variable radius and number concentration. These are (i) nucleation mode (number median radius r<5 nm, only soluble), (ii) Aitken mode (r from 5 to 50 nm), (iii) accumulation mode (r from 50 to 500 nm), and (iv) coarse mode (r>500 nm) (Stier et al., 2005; Zhang et al., 2012). Interactions with clouds are implemented through an explicit activation scheme based on Köhler theory (Abdul-Razzak and Ghan, 2000), with an empirical estimation of maximum supersaturation derived from explicit parcel model calculations. Activated droplet numbers are passed on to a two-moment cloud microphysics scheme (Lohmann et al., 2007; Lohmann and Hoose, 2009) with prognostic variables for cloud droplet number concentration (CDNC) and ice crystal number concentration (ICNC). Emissions and dry and wet deposition are handled consistently between the aerosol scheme and the gas-phase chemistry scheme MOZ (see Sect. 2.3).

HAM can be run either with or without the detailed gas-phase chemistry scheme of MOZ. If run without, then climatological fields from MOZ are used to prescribe monthly mean mixing ratios of oxidants, i.e., ozone, OH, H₂O₂, NO₂, and NO₃. If run interactively, the surface areas of HAM aerosols are used as input for the calculation of heterogeneous reaction rates (Stadtler et al., 2017). At present there is no interaction of aerosols with gas-species photolysis (MOZ uses a climatology of aerosols and lookup tables). These interactions were found to be negligible in an earlier study (Pozzoli et al., 2008a).

2.3 MOZ1.0

The Jülich Atmospheric Mechanism (JAM) 002, which forms the basis of MOZ1.0, has its foundation in a blend of the stratospheric chemistry scheme of the Whole Atmosphere Chemistry Climate Model (WACCM; Kinnison et al., 2007) and the tropospheric Model of Ozone and Related Tracers (MOZART) version 4 (Emmons et al., 2010). The combined chemistry scheme of WACCM and MOZART has been enhanced with a detailed representation of the oxidation of isoprene following the Mainz Isoprene Mechanism 2 (MIM-2; Taraborrelli et al., 2009) and by adding a few primary volatile organic compounds and their oxidation chains. The isoprene oxidation scheme includes recent discoveries of 1,6 H-shift reactions (Peeters et al., 2009), the formation of epoxide (Paulot et al., 2009), and the photolysis of HPALD (Wolfe et al., 2012). Some of the reaction products and rates were taken from the Master Chemical Mechanism version 3.3.1 (Jenkin et al., 2015).

Table 1 lists the primary volatile organic compounds together with their respective oxidants. Radical–radical reactions have been substantially revised since (Emmons et al., 2010). In contrast to the Master Chemical Mechanism, MOZART and JAM002 do not use a radical pool, but instead follow the pathways of peroxy radical reactions with HO₂, CH₃O₂, and CH₃COO₂ (peroxy acetyl) as explicitly as possible. Inorganic tropospheric chemistry considers ozone, NO, NO₂, NO₃, N₂O₅, HONO, HNO₃, HO₂NO₂, HCN, CO, H₂, OH, HO₂, H₂O₂, NH₃, chlorine and bromine species, SO₂, and oxygen atoms.

Six heterogeneous reactions are considered in the troposphere. These are

1. uptake of ozone on dust, including the formation of HO₂;

2. uptake of HO₂ on aqueous aerosol, including cloud droplets and yielding H₂O₂;

3. uptake of NO₃;

4. uptake of NO₂;

5. uptake of HNO₃ on sea salt and dust; and

6. uptake of N₂O₅.

Details on the heterogeneous reactions in ECHAM-HAMMOZ and a discussion of their relevance are given in Stadtler et al. (2017).

The stratospheric chemistry scheme explicitly treats the oxidation and photolysis of 21 halogenated compounds listed in Table 2 together with their approximate lifetimes. Heterogeneous reactions occur on four types of particles:

1. liquid binary sulfate (LBS);

2. supercooled ternary solution (STS);

3. nitric acid tri-hydrate (NAT); = and

4. water ice.

For details, see the Supplement of Kinnison et al. (2007).

The complete listing of the JAM002 species and reaction equations can be found in Supplement 1 (Tables S1–S23). In total, JAM002 contains 246 species and 734 reactions, including 142 photolysis reactions.

MOZ uses the same chemical preprocessor as CAM-Chem (Lamarque et al., 2012) and WACCM (Kinnison et al., 2007) to generate FORTRAN code, which contains the chemical solver for a specific chemical mechanism. In ECHAM-HAMMOZ, all reactions are treated with the semi-implicit (Euler backward integration) solver. This solver uses efficient sparse matrix techniques (LU decomposition and Newton–Raphson iteration) and is set up as follows: within the outer time step loop, up to 11 iterations are performed to achieve a solution within the prescribed relative accuracy. For ozone, NO, NO₂, NO₃, HNO₃, HO₂NO₂, N₂O₅, OH, and HO₂, a relative error of less than 10⁻⁴ is required and less than 10⁻³ for all other species. If convergence is not reached after 11 iterations, the time step is halved and the calculation is repeated. This may happen up to five times. If convergence is still not achieved, a warning message is written into the log file, and the calculation continues. A 3-day test simulation with detailed diagnostics on the solver behavior showed no cases in which the time step length had to be reduced, and convergence was always reached after two to six iterations. As expected, the largest number of iterations occurred under conditions of sunrise and sunset. The model is parallelized in a way that blocks of entire vertical columns on several adjacent longitudes are passed to the solver together, and the convergence threshold is evaluated for the entire block for efficiency reasons. This implies that changing the vector length of the parallelization will affect the results of the chemical calculations (within the error limits given above). More details on the MOZ chemical solver can be found in the Supplement of Kinnison et al. (2007).

The preprocessor code is available with the ECHAM-HAMMOZ model distribution. A simplified chemistry scheme for stratospheric applications (GEOMAR Atmospheric Mechanism; GAM) is also available and can easily be used in lieu of the extensive JAM002 mechanism (in preparation).

Table 1 Primary volatile organic compounds and their oxidants in the JAM002 mechanism. BIGALKANE is a lumped species for all alkanes C₄ and greater, BIGENE lumps all alkenes C₄ and greater, and o-, m-, and p-xylene are lumped into one xylene species. CH₄ is also oxidized by O(¹D) and F, CH₂O is also oxidized by O and HO₂, DMS is also oxidized by BrO, and CH3Br and HCN are also oxidized by O(¹D).

Table 2 Halogenated compounds in JAM002 with relevance for stratospheric ozone chemistry and their approximate lifetimes in years (Miller et al., 1998; Liang et al., 2010; WMO, 2014; Harrison et al., 2016).

3 Chemical parameterizations

3.1 Emissions

All emissions in the ECHAM-HAMMOZ model are controlled via a single “emi_spec” file, which provides a simple and compact way to define all trace gas and aerosol emissions used in a model simulation and ensures proper documentation of the emissions used in a specific run. The emi_spec file consists of three sections: (1) definition of emission sectors and the source of emission information for this sector; (2) the species–sector matrix controlling which emission sectors are active for which species; and (3) an alias table that allows the mapping of species names from emission files to the names that are defined in the chemical mechanism. The emi_spec file that was used in the simulations of this paper is provided in Supplement 2.

In the sector definition, users can specify if emissions from that sector shall be read from the file or if an interactive parameterization (if available) shall be applied. In addition, it is possible to specify a single number to be used as a globally uniform emission mass flux. Furthermore, it can be decided to apply the emissions as a boundary flux condition to the lowest model level, to inject them at the model level near 50 m of altitude (smoke stack emissions), or to distribute them within a specific range of the atmosphere. For biomass burning emissions a special option is available to distribute them across the planetary boundary layer. Finally, the user can also select if emissions shall be interpolated in time or not and if the year of the date information in an emissions file shall be ignored in order to treat emissions as climatology. The simulations shown in this article were performed without time interpolation and using emissions for specific years.

The species–sector matrix has a single float number or a dash in each cell. The float number can be used to easily scale emissions from a particular sector for a particular species or to define emissions for compounds for which no emission data are available by scaling these emissions to those of other compounds (this requires an entry in the alias table). A dash indicates that no emissions for this compound are available in the given sector and is distinct from a zero, which would attempt to read or calculate emissions and then scale them to zero afterwards.

Emission files can be provided in any time resolution (minimum daily). Normally, all emissions files contain monthly data, except for fire emissions, which are provided in daily resolution in the standard configuration.

Table 3 Emissions of trace gases and aerosols used in the standard configuration of ECHAM-HAMMOZ for the year 2008.

With JAM002 and either HAM-M7 or HAM-SALSA as an aerosol module, ECHAM-HAMMOZ has emissions for a total of 43 species that are emitted in 20 sectors. Table 3 lists all emissions for the year 2008.

In the standard configuration of ECHAM-HAMMOZ, the following emissions are calculated interactively:

1. VOC emissions from terrestrial vegetation (MEGAN; Guenther et al., 2012; implementation details in Henrot et al., 2017);

2. DMS emissions from the oceans (Kloster et al., 2006; Lana et al., 2011);

3. dust (Tegen et al., 2002; Stier et al., 2005);

4. sea salt (Guelle et al., 2001); and

5. volcanic sulfur (Dentener et al., 2006b).

Emissions from agriculture (AGR) and waste burning (AWB), forest fires (FFIRE) and grassland fires (GFIRE), aircraft (AIRC), domestic fuel use (DOM), energy generation, including fossil fuel extraction (ENE), industry (IND), ship traffic (SHP), solvent use (SLV), transportation (TRA), and waste management (WST) are taken from the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP; Lamarque et al., 2010). More specifically, as the simulations described here focus on the period after 2000, we make use of the Representative Concentration Pathway (RCP) 8.5 emissions (van Vuuren et al., 2011). The original files, which had a temporal resolution of 10 years, were interpolated to individual years and seasonal cycles were added (Granier et al., 2011). The netCDF emission data files are available from a WebDAV server at the Forschungszentrum Jülich (http://accmip-emis.iek.fz-juelich.de/) and contain detailed README and metadata information.

Ocean emissions of reactive VOCs were obtained from the POET project (Granier et al., 2005), and terrestrial DMS emissions are from Dentener et al. (2006b).

3.2 Lightning

As described in Rast et al. (2014), lightning NO_x emissions are parameterized as a function of the average convective updraft velocity $\stackrel{\mathrm{‾}}{w}$ in a model column following Grewe et al. (2001). Flash frequency is calculated as

$$\begin{array}{}\text{(1)}& F=a(\stackrel{\mathrm{‾}}{w}/{\stackrel{\mathrm{‾}}{w}}_{\mathrm{0}}\sqrt{d/d_{\mathrm{0}}}{)}^{\mathit{\beta }},\end{array}$$

with $a=\mathrm{1.54}\times {\mathrm{10}}^{-\mathrm{5}}$, β=4.9, ${\stackrel{\mathrm{‾}}{w}}_{\mathrm{0}}=\mathrm{1}$ m s⁻¹, and d₀=1 m; d is the vertical cloud extent. Due to the coarse model resolution, a correction factor of 0.7 is applied to the result of this formula to yield a global flash frequency of 49 flashes per second for the year 2008, which is within the uncertainty of 44 ± 5 flashes per second observed from the Optical Transient Detector satellite instrument during 1995 to 2000 (Christian et al., 2003).

The fraction of cloud-to-ground flashes is calculated according to Price and Rind (1994) as

$$\begin{array}{}\text{(2)}& f_{\mathrm{cg}}=(\mathrm{0.021}{d}_{\mathrm{c}}^{\mathrm{4}}-\mathrm{0.648}{d}_{\mathrm{c}}^{\mathrm{3}}+\mathrm{7.49}{d}_{\mathrm{c}}^{\mathrm{2}}-\mathrm{36.54}{d}_{\mathrm{c}}{)}^{-\mathrm{1}},\end{array}$$

where d_c denotes the cold cloud thickness, i.e., the vertical extent of the part of the cloud with temperatures below freezing. Following Price et al. (1997) the amount of NO generated per flash is given as 1×10¹⁷ $\mathrm{molec}.{\mathrm{J}}^{-\mathrm{1}}$, and average flash energies are assumed to be 6.7×10⁹ J for cloud-to-ground flashes; one-tenth of this is for intra- and inter-cloud flashes. With these factors applied, the global amount of NO generated from lightning in the year 2008 would be 5.05 Tg(N), which is well within the range of other estimates (e.g., Schumann and Huntrieser, 2007) and was recommended as a target value in earlier model intercomparison projects. As we found a significant influence of global lightning NO_x on global tropospheric ozone and OH (see Sect. 6) with methane and methylchloroform lifetimes more in the range of the literature values at lower lightning NO_x, we scaled the lightning emissions down. In the default configuration of the model, global lightning NO_x emissions from the simulation described in this study for the year 2008 are 1.2 Tg(N). Within the model column, the NO from lightning is distributed according to the climatological vertical profiles of Pickering et al. (1998).

3.3 Lower boundary conditions for long-lived stratospheric species

Halogenated species, which are primarily decomposed in the stratosphere, are not emitted into the model, but instead a lower boundary condition is specified for these compounds. In addition to the species in Table 2 the model also specifies lower boundary conditions for N₂O, CH₄, and CO₂. The latter can be turned off if the model is run with all carbon cycle components. The lower boundary conditions are provided as zonally averaged, monthly values from the Whole Atmosphere Chemistry Climate Model (WACCM) input for the simulations in the Chemistry–Climate Model Initiative (CCMI) initiative (see Sect. 2.3.2 in Tilmes et al., 2016). The organic halogen scenario (here, RCP8.5) is based on WMO (2011) and described in Eyring et al. (2013) and Morgenstern et al. (2017). Boundary conditions for N₂O, CH₄, and CO₂ are taken from Meinshausen et al. (2011) as recommended by CCMI. As described in Tilmes et al. (2016), the boundary conditions used in ECHAM-HAMMOZ include a latitudinal gradient based on aircraft measurements from the HIAPER (High-Performance Instrumented Airborne Platform for Environmental Research) Pole-to-Pole Observations (HIPPO) campaigns (Wofsy et al., 2011).

The lower boundary condition for methane uses a relaxation to the climatological values with an e-folding time of 3 days in order to preserve regional methane emission patterns while at the same time preventing a drift of methane concentrations due to possibly unbalanced sources and sinks (see Rast et al., 2014).

3.4 Photolysis

As described in the Supplement of Kinnison et al. (2007), photolysis frequencies are calculated as a product of the prescribed extra-atmospheric solar flux, the atmospheric transmission function (dependent on model-calculated ozone and O₂), the molecular absorption cross section, and the quantum yield of the specific reaction channel. There are 34 channels in the wavelength regime from 120 to 200 nm and 122 channels between 200 and 750 nm. At wavelengths less than 200 nm, the transmission function is computed explicitly, and absorption cross sections and quantum yields are prescribed, except for O₂ and NO, for which detailed parameterizations are used (see the Supplement of Kinnison et al., 2007). Beyond 200 nm, the transmission function is calculated from a lookup table as a function of altitude, column ozone, surface albedo, and zenith angle. The maximum zenith angle for which photolysis frequencies are calculated is 97^∘. The temperature and pressure dependence of absorption cross sections and quantum yields is also interpolated from lookup tables.

The UV albedo is parameterized according to Laepple et al. (2005) using satellite-derived albedo maps for snow and non-snow conditions based on the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument. The threshold for switching from non-snow to snow values is a snow depth of 1 cm calculated by ECHAM. Sea ice albedo is prescribed as 0.78 in the Northern Hemisphere and 0.89 in the Southern Hemisphere. The albedo over ice-free oceans is set to 0.07. The influence of clouds is parameterized according to Brasseur et al. (1998) through the computation of an effective albedo and modification of the atmospheric transmission function. Both effects are combined into a single factor that varies by model level.

3.5 Dry deposition

Deposition of trace gases and aerosols on the Earth's surface is parameterized according to the resistance model of Wesely (1989) using a big-leaf approach for vegetated surfaces (Ganzeveld and Lelieveld, 1995). The ECHAM-HAMMOZ dry deposition scheme distinguishes between vegetated land surfaces, bare soils, water, and snow or ice and uses the corresponding surface types from the JSBACH land model (see Sect. 2.1).

As described in Stier et al. (2005), the leaf area index is taken from an NDVI (normalized differential vegetation index) climatology (Gutman et al., 1995) and the Olson (1992) ecosystem database. Canopy height, roughness length, and forest fraction are from a climatology. Soil pH is specified for 11 different soil types. Technically, these parameters can also be obtained online from the JSBACH land surface model. This has been tested in Stanelle et al. (2014), but this coupling has not been used in the simulations described in this study.

Surface resistances are explicitly specified for H₂SO₄, HNO₃, NO, NO₂, O₃, and SO₂. Those of other species are calculated relative to O₃ and SO₂ following Wesely (1989). A dynamic coupling of stomatal resistance with JSBACH has not been implemented. Henry coefficients and reactivity factors have been defined for 137 species (see Table S24 in Supplement 1). All of these species can be deposited; however, for most of them the deposition rates will be very small due to low Henry values or low reactivities. Where available, Henry coefficients were taken from the literature, and in other cases we assumed that molecules with similar structures have Henry values on the same order of magnitude. In particular, OOH groups were considered similar to OH groups in terms of their impact on Henry constants. Organic peroxides (ROOH) are assumed to have a surface reactivity f₀ of 1. Other organic molecules with oxygen were assigned with f₀=0.1. Note that while HO₂ can undergo dry deposition (H=690 M atm⁻¹), we did not define Henry values and surface reactivities for organic peroxy radicals, although some of these (notably from higher oxidation states) might well be soluble and could therefore be efficiently removed by dry deposition.

3.6 Wet deposition and scavenging

The ECHAM-HAMMOZ wet deposition scheme is based on Croft et al. (2009) for below-cloud scavenging by rain and snow and Croft et al. (2010) for in-cloud (nucleation and impaction) scavenging. The in-cloud scavenging scheme treats nucleation and impact scavenging in stratiform and convective clouds and distinguishes between warm, cold, and mixed-phase clouds. For aerosols, scavenging also takes place below clouds in rain and snow. For gases, the gas-phase liquid-phase equilibrium is calculated based on Henry's law (see Table S2), and no below-cloud scavenging takes place except for HNO₃ and H₂SO₄.

4 Simulation setup

The simulations described in this article are based on the ECHAM-HAMMOZ model version and input datasets as released in February 2017. The base run is a 10-year simulation from October 2002 through December 2012. The simulation started with initial conditions for October 2002 taken from the Monitoring Atmospheric Composition and Climate (MACC) reanalysis described by Inness et al. (2013). The first 3 months are used as spin-up time and are not used in the analysis of results. Note that the stratosphere takes about 4 years to reach a dynamically stable state (see Fig. S3.1 in the Supplement 3).

The model resolution is $\mathrm{1.875}{}^{\circ }\times \mathrm{1.875}{}^{\circ }$ in longitude and latitude (spectral triangular truncation T63), with 47 vertical levels from the surface to 0.01 hPa (full level pressure). The corresponding model time step is 7.5 min. Sea surface temperatures and sea ice coverage are prescribed for each year of simulation, following the Coupled Model Intercomparison Project Phase 5 (CMIP5) AMIP simulation protocol (Giorgetta et al., 2012). In addition, temperature, vorticity, divergence, and surface pressure from 6-hourly output of the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis (Dee et al., 2011) were used to specify the dynamics of “real weather” in nudging mode (Kaas et al., 2000). The simulation includes full interaction between aerosols and gas-phase chemistry, full feedback of aerosols and trace gases on the radiation, and feedback of chemically produced water in the stratosphere onto the climate model. All simulations performed for this study use the HAM-M7 aerosol scheme for aerosol formation and microphysical processes.

Some of the analyses presented below discuss the 10-year core period of the simulation, while the more detailed comparisons with observations focus on the year 2008. This year was chosen as a reference year by the HAMMOZ consortium because of data availability and because it has been used in several other studies and in the Copernicus Atmospheric Monitoring Service validation activities (Eskes et al., 2015). The year 2008 began with La Niña conditions and turned neutral by June. The global average temperature was 0.49 K above the 20th century mean (NOAA, 2008), but looking back from 2012, the end year of our simulation, it is not among the 10 warmest years on record (NOAA, 2012). In terms of Antarctic ozone depletion, 2008 was relatively close to the 1979–2007 average (WMO, 2008).

Table 4 Summary of simulations performed for the analysis of global trace gas budgets in this article.

As described above, the released model version has very low lightning NO emissions, and the parameterization of the heterogeneous reaction of HNO₃ neglects the potential of reevaporation of aerosol nitrate to gaseous HNO₃ (Stadtler et al., 2017), which has an impact on the total amount of deposited nitrogen and also affects the budgets of ozone and nitrogen oxides as shown below. In order to investigate the impacts of these two factors we performed a small series of sensitivity simulations for the year 2008 based on restart files of the reference run. Table 4 briefly describes these simulations. Two of the simulations double and quadruple the amount of lightning emissions, respectively, so that they are more in line with current estimates ranging from 2 to 8 Tg(N) yr⁻¹ (Schumann and Huntrieser, 2007). The other simulation deactivates the heterogeneous reaction of HNO₃ on sea salt and dust and is otherwise identical to the base run.

5 Comparison with observations

5.1 Total column ozone and stratospheric processes

We begin our model evaluation with a discussion of seasonal total column ozone (TCO) distributions in comparison with retrievals from the Infrared Atmospheric Sounding Interferometer (IASI), onboard the sun-synchronous polar-orbiting MetOp platforms, and from the Ozone Monitoring Instrument (OMI) onboard the Aura satellite (Levelt et al., 2006). Furthermore, the 10-year evolution of TCO has been compared to the The Modern-Era Retrospective Analysis for Research and Applications version 2 (MERRA-2), for which a large variety of satellite retrievals were assimilated (Gelaro et al., 2017).

IASI is a high-resolution Fourier transform spectrometer designed to measure the outgoing thermal infrared radiation from the Earth's surface and the atmosphere in a nadir-viewing geometry (Clerbaux et al., 2009). The first IASI instrument was launched in 2006 on the MetOp-A platform and it is still operating. A second instrument onboard MetOp-B was launched in 2012. IASI observations – a set of four simultaneous footprints of 12 km at nadir – are performed every 50 km along the track of the satellite at nadir and across-track over a swath width of 2200 km. IASI crosses the Equator at around 09:30 and 21:30 mean local solar time, achieving near global coverage twice a day. The vertical abundance and distribution of O₃ and CO are retrieved in near real time from individual IASI measurements with the Fast Optimal Retrievals on Layers for IASI (FORLI) algorithm (fully described in Hurtmans et al., 2012). IASI measurements are filtered out based on the cloud coverage and several quality flags to ensure stable and consistent products. The estimated errors on the retrieved total columns are mainly latitude dependent. For O₃, such errors are usually below 3 %, but slightly larger (∼ 7.5 %) when the signal is particularly weak over the Antarctic ozone hole or due to strong water vapor influence at tropical latitudes (Hurtmans et al., 2012). The error on the retrieved CO total columns is generally below 10–15 % at middle and tropical latitudes, but increases up to 20–25 % in polar regions (George et al., 2015). All details as to the retrieval methodology, characterization of the retrieved products, and validation against a large suite of independent ground-based, airborne, and satellite measurements can be found in previous studies (Hurtmans et al., 2012; George et al., 2015; Boynard et al., 2016, 2018; Wespes et al., 2016, 2017).

For a meaningful model-to-satellite comparison, daily average ECHAM-HAMMOZ vertical profiles of O₃ and CO have been smoothed with the IASI averaging kernels (A) to account for the vertical sensitivity of the instrument to the species. The model profiles were firstly vertically interpolated onto the IASI retrieval levels. Then the smoothing of the model profiles (X_model) to the lower vertical resolution of the IASI measurements was performed following the formalism of Rodgers (2000):

$$\begin{array}{}\text{(3)}& X_{\mathrm{smoothed}}=X_{\mathrm{a}}+A\cdot (X_{\mathrm{model}}-X_{\mathrm{a}}),\end{array}$$

where X_smoothed is the smoothed model profile and X_a is the a priori profile. As numerous IASI observations are usually available for a given day and grid cell, the total columns calculated from the smoothed model profiles were averaged to generate one daily mean value.

Figure 1 compares seasonal mean total ozone columns of the year 2008 from IASI with the ECHAM-HAMMOZ outputs smoothed by the IASI sensitivity functions. The model generally captures the highs and lows that are observed by the satellite as well as the changes in these patterns throughout the seasons. It shows a tendency to underestimate the total ozone column in middle to high northern latitudes (< 10 %; i.e., 5–30 DU) and over the Antarctic (< 15 %; i.e., 10–40 DU) in all seasons, while it overestimates over the tropical regions (< 5 %; i.e., 10 DU) and the southern midlatitudes (∼ 10 %; i.e., 25 DU). These IASI–model discrepancies are generally larger than the errors on the retrieved TCO (see above) and the IASI-UV instrument biases recently reported in the validation experiment of Boynard et al. (2018). These authors report errors that are generally lower than 5 %, with a global bias of 0.1–2 %.

Figure 1 (a, b) Comparison of seasonal mean total column ozone (in DU) from ECHAM-HAMMOZ with data from the IASI instrument. Data and model results are from 2008. IASI data were interpolated onto the model grid and averaging kernels were applied to the model data. (c) TCO relative differences (in %) between ECHAM-HAMMOZ and IASI.

Figure 2 Comparison of the seasonal cycle of zonal mean total column ozone (in DU) between ECHAM-HAMMOZ (a) and Aura OMI (b). Data and model are for the year 2008. Aura OMI data were interpolated to model grid.

OMI is an ultraviolet visible nadir-viewing solar backscatter spectrometer. We use the level 3 data product, which is globally gridded to 1^∘ latitude by 1^∘ longitude. The Aura OMI data span the temporal range from 2004 to present. Daily OMI data were interpolated to the model resolution and compared to the model results as latitude–time series plots in Fig. 2. Unlike the comparison of the model to IASI observations in Fig. 1, this analysis has not applied the OMI averaging kernel to the model results, and therefore one needs to be careful in overstating any biases between the model and the observations. A discussion of biases among satellite TCO products is beyond the scope of this paper.

The purpose of Fig. 2 is to document that the model nicely represents the overall zonal mean latitude versus time structure compared to OMI. Both the Northern and Southern Hemisphere spring maxima are consistent with OMI along with the tropical minimum in TCO around 260 DU. The Southern Hemisphere ozone hole is also well represented, suggesting that the heterogeneous chemistry defining the Antarctic ozone hole is adequately formulated in the model (see Solomon, 1999, and references within).

Figure 3 The 10-year latitude–time variations in the total ozone column in ECHAM-HAMMOZ (a) and in the MERRA-2 reanalysis (b). A difference plot can be found in the Supplement.

Figure 3 shows a comparison between the latitude–time variations in TCO during the 10 years of our ECHAM-HAMMOZ reference simulations with TCO from the MERRA-2 reanalysis (Gelaro et al., 2017). Generally, the differences in TCO remain below 6 % (see also Fig. S3.1 in the Supplement), which is similar to the accuracy estimate of Wargan et al. (2017) for the MERRA-2 stratospheric ozone field. The interannual variabilities of the springtime ozone maximum in the Northern Hemisphere and the Antarctic minimum are well captured. The model slightly underestimates TCO in the tropics, and it has a tendency to underestimate the Southern Hemisphere midlatitude maxima during austral spring. Closer inspection of the differences (Fig. S3.1) reveals larger biases during the first 4 years of the simulation. These can be attributed to the initial conditions, which, although suitable for October 2002, i.e., the beginning of our simulation, were taken from a different model and might therefore reflect a somewhat different stratospheric state than what ECHAM-HAMMOZ would have simulated with a longer spin-up period. A decay time of 4 years for TCO disturbances is quite reasonable and indicates a reasonable representation of the stratospheric age of air, which has not been explicitly evaluated for this model.

Figure 4 Comparison of temperature, HNO₃, and O₃ between Aura MLS observations (a) and ECHAM-HAMMOZ (b). Data and model are for the year 2008.

Figure 4 further explores the representation of the Antarctic region by showing a time-dependent vertical cross section of key constituents at 81^∘ S. In this figure, model results of temperature, nitric acid (HNO₃), and ozone are compared to daily binned (4.5^∘ latitude × 10^∘ longitude) data from version 4 of the Microwave Limb Sounder (MLS) instrument onboard the Aura satellite. Details of the accuracy and precision of the MLS observations are discussed in Livesey et al. (2016). Since polar heterogeneous chemistry is very temperature sensitive (e.g., Solomon et al., 2015) it is important to have an accurate representation of temperature when comparing model results to observational data for a given year. Even in a nudged simulation, temperature biases in the lower stratosphere can amount to a few degrees (Krefting, 2017). Figure 4, row 1 shows excellent agreement between the retrieved MLS temperatures and those used in ECHAM-HAMMOZ, giving confidence that there are no model temperature biases that could impact the heterogeneous chemistry derived in the model.

Another important constituent to show is HNO₃. Here, the HNO₃ gas-phase abundance is affected by the formation of NAT PSC particles, which can settle out of the stratosphere and cause irreversible denitrification. Figure 4, row 2 shows that the model adequately represents the HNO₃ abundance and the process that controls the loss of total inorganic nitrogen in the model atmosphere. If anything, the model over-denitrifies by 0.5 nmol mol⁻¹. This result is consistent with use of the current ECHAM-HAMMOZ heterogeneous chemistry module (Kinnison et al., 2007) in other model frameworks, e.g., the Whole Atmosphere Chemistry Climate Model version 3 (see evaluation in SPARC, 2010). The model does show a low bias in HNO₃ between January and May in the pressure range of 80–10 hPa. This low bias is not due to the effect of denitrification, since this process starts in June. The model does not include a representation of lower mesospheric and stratospheric particle production of NO_x. However, the descent of NO_x from these processes and the eventual formation of HNO₃ would not affect the lower stratosphere until later in the austral winter. The model does include the heterogeneous conversion of N₂O₅ to HNO₃ on sulfate aerosol (Austin et al., 1986). Future work to better understand this bias will need to focus on the tropical gas-phase net production of NO_y and subsequent meridional transport to higher latitudes.

The ozone evolution is shown in Fig. 4, row 3. The overall representation of ozone from the lower mesosphere to the lower stratosphere is captured by the model. However, in the lower stratosphere spring, the model is biased high by approximately +0.5 µmol mol⁻¹. This high ozone bias translates to the model having too much TCO in the Antarctic spring period, which is consistent with Fig. 2.

5.2 Tropospheric ozone

Figure 5 Annual mean bias of ECHAM-HAMMOZ with respect to ozonesonde data from Tilmes et al. (2012) at 250 hPa (a), 500 hPa (b), and 800 hPa (c). At 250 hPa the relative error is shown for better visibility. Figure layout and scales are comparable to Lamarque et al. (2012).

Figure 6 Frequency distributions of the ozone bias (in nmol nmol⁻¹) at the 42 stations from Tilmes et al. (2012). Note: high northern latitude stations with biases larger than ±40 nmol nmol⁻¹ at 250 hPa are not shown.

Figure 5 shows the year 2008 annual mean bias of ozone in comparison to the ozonesonde climatology of Tilmes et al. (2012) at three different pressure levels. For ease of comparison, we have chosen a similar layout and scale as Lamarque et al. (2012). At 250 hPa, the biases are generally between −35 and +35 nmol nmol⁻¹ with the exception of high northern latitude stations, where the bias is as low as −114 nmol nmol⁻¹ (Eureka and Resolute, Canada) and Prague, Czech Republic, where the bias is +66 nmol nmol⁻¹. The model overestimate at high northern latitudes is qualitatively similar to CAM-Chem (Lamarque et al., 2012). We concur with Lamarque et al. (2012) that the reason is likely associated with a mismatch between the model tropopause and the real tropopause in this region. Due to the very steep gradients of ozone around the tropopause, even small vertical displacements can lead to large discrepancies in simulated ozone values if the comparison is made on pressure levels. Future work should probably consider evaluating models with ozonesonde data relative to the tropopause. For comparison, Fig. S3.2 in the Supplement presents plots similar to Fig. 5, but for the simulation lght*4.

In contrast to CAM-Chem, the Northern Hemisphere midlatitude biases in ECHAM-HAMMOZ are more or less equally distributed. One may discern a tendency of the model to overestimate ozone at 250 hPa around the Pacific, whereas there appears to be underestimation around the Atlantic.

At 500 hPa, the pattern of the ECHAM-HAMMOZ biases is similar to that at the 250 hPa level, but the values are generally smaller. Only 6 stations out of 42 have biases with absolute values larger than 10 nmol nmol⁻¹, and 5 of these 6 stations are located in the tropics. Ascension, Natal, and Reunion exhibit large low biases, whereas high biases are found at Hilo, Hawaii, and Samoa. At 800 hPa the biases are somewhat shifted to more positive values so that no site has a low bias of more than 10 nmol nmol⁻¹, and four stations show a high bias greater 10 nmol nmol⁻¹. Averaged over all 42 sites, the mean biases at 250, 500, and 800 hPa are −17.8, −1.5, and +1.6 nmol nmol⁻¹, respectively. If we exclude the high-latitude Northern Hemisphere stations, the bias at 250 hPa is reduced to −3.6 nmol nmol⁻¹. Figure 6 shows frequency distributions of the model biases at the three pressure levels.

The seasonal cycle of tropospheric ozone is evaluated with the help of Taylor plots in Fig. 7. Similar to Lamarque et al. (2012) we show regional averages at the three model levels of Fig. 5. However, we retain the original region definitions and color codings of Tilmes et al. (2012). Taylor plots with individual stations and also for the sensitivity run lght*4 can be found in Fig. S3.3. For seven out of nine regions, the correlation between the observed and simulated seasonal cycle is positive so that the symbols appear in Figure 7. Exceptions are Canada and the tropics. The normalized root mean square error (concentric gray circles around a standard deviation ratio of 1 and a correlation of 1) is generally below 0.8. Exceptions are the eastern Northern Hemisphere polar stations at 250 hPa, the Southern Hemisphere polar stations at 500 and 800 hPa, and the Southern Hemisphere midlatitude stations at 800 hPa. At 250 hPa the correlation is better than 0.7 over most regions, and exceptionally good results are obtained over the Northern Hemisphere western polar region. The correlation slightly worsens at 500 and 800 hPa, but generally remains better than 0.6. Across all 42 sites, the average correlation coefficients at 250, 500, and 800 hPa are 0.59, 0.59, and 0.68, respectively. If we leave out the tropics, which have the worst correlation, they increase to 0.68, 0.68, and 0.73, respectively.

Hence, as a summary, we can state that tropospheric ozone in the base run is very well simulated with two exceptions: (1) there is a severe underestimate at high northern latitude sites at 250 hPa, and (2) the (small) seasonal cycle over the tropics is not well captured.

Figure 7 Taylor plots of regional averaged ozonesondes at 250 hPa (a), 500 hPa (b), and 800 hPa (c). Ozonesonde data and region definitions are from Tilmes et al. (2012). Where a region is not shown in a panel, the respective data point is outside the axis range.

5.3 Surface ozone

Figure 8 Comparison of gridded rural surface ozone observations from the TOAR database (Schultz et al., 2017) (a) with ECHAM-HAMMOZ output (b) for January and July 2008. The model results have been regridded to 5^∘ × 5^∘ to match the resolution of the observations.

Figure 9 Monthly mean bias of surface ozone mixing ratios in January and July 2008 for all 5^∘ × 5^∘ grid boxes for which the TOAR database contains data in 2008.

For the evaluation of ECHAM-HAMMOZ with surface ozone observations, we make use of the recently published database from the Tropospheric Ozone Assessment Report (TOAR). As described in Schultz et al. (2017), this database contains hourly data from more than 9000 scientific and air quality monitoring stations worldwide, and it has a globally consistent classification scheme to distinguish urban from rural locations. The classification scheme is based on threshold combinations of global satellite data products of nighttime light intensity, population density, and OMI NO₂ columns. For details see Schultz et al. (2017).

Figure 8 shows gridded maps of TOAR data at rural stations (top row) in comparison with ECHAM-HAMMOZ base run output regridded to the same resolution of 5^∘ × 5^∘ for January and July 2008. Figure 9 displays the bias between the model and observations. The first thing that becomes apparent in Fig. 8 is the evident geographical bias of the observation database. About three-quarters of the grid boxes with measurements by rural stations are located either in Europe or North America, and the rest is scattered across the world. The problem of insufficient observational coverage of reactive gas measurements is widely known, and the community has yet to develop a sound strategy to deal with it.

Where measurements exist, the model generally shows good agreement with the observations in both January and July. During the boreal summer, ozone over the eastern US and the North Sea–Baltic Sea region is somewhat overestimated. Closer inspection reveals differences of up to −25 and +30 nmol nmol⁻¹ in individual grid boxes, for example over the Mediterranean or in Nebraska, US. However, altogether the model yields excellent agreement with mean bias of 1.13 nmol mol⁻¹ in January and 5.28 nmol mol⁻¹ in July. Additional information, also concerning the sensitivity experiments lght*2 and lght*4, can be found in Figs. S3.4 and S3.5.

5.4 Total column CO

Figure 10 (a, b) Comparison of seasonal mean total column CO (in 10¹⁸ cm⁻²) by ECHAM-HAMMOZ with data from the IASI instrument. Data and model results are from 2008. IASI data were interpolated onto the model grid and averaging kernels were applied to the model data. (c) CO total column relative differences (in %) between ECHAM-HAMMOZ and IASI.

Figure 10 shows seasonal mean total column CO from IASI in comparison with ECHAM-HAMMOZ smoothed by the IASI averaging kernels. The model reproduces the seasonal variations with the main features of the retrieval, but also shows a couple of differences.

1. The model tends to underestimate CO by about 15 % in the Northern Hemisphere, especially during spring and summer. Apparently, ECHAM-HAMMOZ loses CO too quickly during these seasons because of relatively high OH levels (see Sect. 6). Inaccurate emissions (see Sect. 5.5) and transport from polluted regions likely also contribute to the differences.

2. In the Southern Hemisphere a positive bias of 15–40 % is found, likely coming from an overestimation of the fire emissions in ACCMIP. The positive bias is particularly pronounced over Indonesia, central Africa, and Amazonia.

5.5 Surface CO

Figure 11 Comparison of monthly mean surface CO measurements from GAW with ECHAM-HAMMOZ base run results for January, April, July, and October 2008. Each symbol represents data from one measurement location.

Figure 11 displays the latitudinal gradients of surface CO concentrations from the World Meteorological Organisation Global Atmosphere Watch (GAW) network (see Schultz et al., 2015) in comparison with ECHAM-HAMMOZ base run results for the year 2008. In general, the model captures the latitudinal variations of CO well throughout the year. However, in higher northern latitudes the simulated CO is underestimated by up to 40 nmol nmol⁻¹ (33 %) in April and to a lesser degree also in January. The reasons for such model–observation discrepancies have been discussed in Stein et al. (2014) and are likely related to inaccurate emissions data in combination with excessive dry deposition of CO. The tendency of ECHAM-HAMMOZ to generate too much OH (see Sect. 6) could also play a role here. The model also overestimates CO in the Southern Hemisphere. This bias is largest during austral winter. The reasons for this bias are unclear at present, but could be related to excessive emissions from biomass burning (see Sect. 5.4). A similar pattern of difference was seen across the suite of ACCMIP models, as described in Naik et al. (2013).

6 Global budgets

6.1 Ozone and OH

Stevenson et al. (2006)Young et al. (2013)Naik et al. (2013)Lamarque et al. (2012)Jöckel et al. (2016)

Table 5 Global tropospheric ozone budgets and tropospheric methane lifetimes of the four ECHAM-HAMMOZ simulations described in Table 4. An ozone threshold of 150 nmol mol⁻¹ is used to denote the tropopause. For comparison, multi-model mean values from Stevenson et al. (2006), Young et al. (2013), and Naik et al. (2013) and the “GEOS5” budget terms from Lamarque et al. (2012) are also included.

^a Computed as whole atmosphere burden of CH₄ over tropospheric loss of CH₄ as in Naik et al. (2013). Values in parentheses were computed according to Jöckel et al. (2016) as tropospheric CH₄ burden over tropospheric CH₄ loss. ^b Values from selected models with relatively low O₃ bias and CH₄ lifetime close to the multi-model mean ^c year 2000 results.^d Burden from 15 ACCMIP models, budget terms from 6 models. ^e Tropopause threshold at O₃<100 nmol mol⁻¹.

The global budgets of ozone and the tropospheric methane lifetime of the ECHAM-HAMMOZ base run are in the range of estimates from other recent models and model intercomparison studies (Table 5). The base run ozone burden is slightly lower than the averages of the multi-model studies. Production and loss terms are above the average plus 1 standard deviation of the selected models in Stevenson et al. (2006), but well in the range of models reported by Young et al. (2013). Since the HAMMOZ chemical mechanism resembles the CAM-Chem mechanism to a substantial degree, one might expect a better agreement with Lamarque et al. (2012), who report substantially lower values. However, Lamarque et al. (2012) used an ozone threshold of 100 nmol mol⁻¹ for the tropopause definition, whereas all other studies in Table 5 used a threshold of 150 nmol mol⁻¹. If we evaluate the base run ozone budget with a 100 nmol mol⁻¹ threshold, we obtain a burden of 292 Tg and production and loss rates of 5192 and 4807 Tg yr⁻¹, respectively.

The above-average ozone production and loss rates are most likely due to the more detailed VOC mechanism in ECHAM-HAMMOZ. Stevenson et al. (2006) noted that earlier model simulations with fewer primary VOC species tended to yield lower production and loss rates. Considering that our base run has very low lightning NO_x emissions (at the low end of the models described in Young et al., 2013), it is somewhat surprising that our ozone chemistry is so active. The ozone lifetime is slightly shorter than the multi-model averages of Stevenson et al. (2006) and Young et al. (2013), but well within the range of these studies. Yet, the ozone lifetime of ECHAM-HAMMOZ is substantially shorter than the result reported by Lamarque et al. (2012) (19.2 days versus 26 days if we use their tropopause threshold of O₃<100 nmol mol⁻¹). As noted by Young et al. (2013): “Despite general agreement on how the drivers impact global-scale shifts in tropospheric ozone, magnitudes of the regional changes and the overall ozone budget vary considerably between different models.”

Figure 12 demonstrates that the tropospheric ozone budget of ECHAM-HAMMOZ remains rather stable over the 10 years of the reference simulation, although some effect from the spin-up can still be seen during the first half of the year 2003. The chemical production and loss terms peak in the boreal summer due to the greater impact of NO_x and VOC emissions from the Northern Hemisphere land masses. We suggest that it might be illustrative to compare the seasonal cycles of individual ozone budget terms from different models in future studies to better understand the reported differences in global annual ozone budgets.

In ECHAM-HAMMOZ the tropical upper troposphere appears to play a prominent role for the global tropospheric trace gas budgets as evidenced by our lightning NO_x sensitivity simulations. Table 6 lists the percent changes in the ozone budget terms when we double or quadruple the lightning NO_x emissions. The tropics are the region with the highest frequency of thunderstorms and the highest flash density (Boccippio et al., 2000), and we do find the largest changes in the ozone budget in this region. For example, the chemical production of ozone increases by 33 % if we increase the lightning NO_x emissions to 4.8 Tg(N) yr⁻¹, a value close to the mean or median lightning NO_x emissions of the models described in Young et al. (2013). The global increase in the ozone production is 1047 Tg yr⁻¹ (Table 5), and the increase in the tropics constitutes 80 % of this change. The Northern Hemisphere is least affected by the increased lightning NO_x source due to the much larger surface and aircraft emissions in this region. This is also reflected in the evaluation of the model runs with surface ozone observations (see Sect. 5.3): in spite of the large changes in the global budget terms and the substantial increase in the global burden, the mean bias in comparison with the gridded TOAR dataset of rural stations increases only moderately from 5.3 to 6.9 nmol mol⁻¹. The density of stations in the Northern Hemisphere is much greater than elsewhere so that the larger changes in surface ozone in the tropics and Southern Hemisphere (see Fig. S3.3) are not accounted for in the bias calculation (Fig. S3.4).

The tropospheric average OH concentration and methane lifetime of our base run are close to the multi-model average diagnosed by Naik et al. (2013). If we increase lightning NO_x emissions, OH increases and the CH₄ lifetime decreases as expected. With 4.8 Tg(N) yr⁻¹ as a global lightning source, the CH₄ lifetime is 8.28 years, which is more than 2 standard deviations below the observational constraints from Prinn et al. (2005) and Prather et al. (2012). However, our CH₄ lifetime appears rather consistent with the lifetime of Jöckel et al. (2016), who use a different method for the calculation (see footnote a of Table 5). Their year 2008 lifetime of 7.65 years falls between our lght*2 and lght*4 simulations if we apply the same method as Jöckel et al. (2016). Given that their lightning NO_x emissions were about 4 Tg(N) yr⁻¹ during this period, the agreement is remarkable. As a consequence we note that chemistry models that use the dynamical core and physics of ECHAM have a tendency to be too reactive and generate too much OH. This has already been an issue in earlier ECHAM-HAMMOZ model versions (e.g., Rast et al., 2014). Indeed, Baumgaertner et al. (2016) have shown that the dynamical core can have a large impact on the global CH₄ lifetime. If we put this information in context with the analysis of regional ozone budget changes due to lightning NO_x emission changes, we hypothesize that there is some issue with the dynamics or physics of ECHAM in the tropical troposphere that impacts its ability to reproduce the global budgets of reactive trace gases. A further hint in this direction is given by Stevens et al. (2013), who pointed out that both ECHAM5 (which forms the basis of the EMAC model reported by Jöckel et al., 2016) and ECHAM6 (the basis of ECHAM-HAMMOZ as described here) have a tropical precipitation bias of up to 5 mm day⁻¹. Figure 13 shows the difference between a decadal ECHAM6.3 simulation (i.e., the ECHAM-HAMMOZ model without the chemistry and aerosol schemes) versus the ERA-Interim reanalysis (Dee et al., 2011), which confirms the statement of Stevens et al. (2013). In addition, the amounts of cloud liquid water and cloud ice are underestimated in the tropics (Lohmann and Neubauer, 2018), pointing to problems either with the parameterization of convection or detrained condensate with implications for the wet scavenging of trace gases and aerosol particles. A more detailed investigation of this issue is beyond the scope of this paper.

Table 6 Percent changes in ozone budget terms in the Northern Hemisphere, the Southern Hemisphere, and the tropics of the simulations lght*2 and lght*4 versus the base run. The latitude boundary of the tropics was chosen at 20^∘ N or S.

Figure 12 Decadal variation in the major tropospheric ozone chemical production and loss terms from the ECHAM-HAMMOZ reference simulation. Blueish colors represent the global ozone production rates (in Tg month⁻¹) due to the reaction of HO₂, CH₃O₂, and other organic peroxy radicals (RO₂); reddish colors show ozone loss rates due to the photolysis and subsequent reaction of O(¹D) with water, due to the reaction of ozone with HO₂, and due to reaction with OH.

Figure 13 Decadal mean bias of total precipitation in ECHAM6.3 compared to the ECMWF ERA-Interim reanalysis (Dee et al., 2011). Gray shading denotes areas that are not significant according to a t test (95 % confidence interval). Figure regenerated from data described in Krefting (2017).

6.2 NO_x budget

The global NO_x budget of the base run differs somewhat from other estimates (e.g., Xu and Penner, 2012) as can be seen from Table 7. While the total NO_x emissions (except for lightning as discussed above) are very similar to other recent studies, the dry and wet deposition rates are about a factor of 2 lower. This is due to the parameterization of the heterogeneous uptake of HNO₃ on sea salt and dust aerosols, which removes HNO₃ from the system and does not allow for reevaporation from the aerosol or droplet phase. The sensitivity run no_het_HNO3 yields dry and wet deposition rates of nitrogen that are very close to Xu and Penner (2012) and other model studies.

Xu and Penner (2012)

Table 7 Global tropospheric NO_x budgets of the four ECHAM-HAMMOZ simulations described in Table 4.

 Xu and Penner (2012) distinguish between gas-phase and aerosol nitrate, while ECHAM-HAMMOZ does not have aerosol nitrate as a separate tracer. The deposition rates listed in this table are total nitrate deposition.

6.3 Radiation, clouds, and aerosol

Even though the focus of this paper is the tropospheric and stratospheric gas-phase chemistry in ECHAM-HAMMOZ, it is useful to evaluate aerosol burdens and cloud and radiation fields in order to assess how the comprehensive gas-phase chemistry mechanism of MOZ relates to HAM-only simulations and other studies. This section focuses on 10-year global mean values of radiation, cloud, and aerosol variables. Spatial distribution maps of speciated aerosol mass burdens and time series of seasonal aerosol burden means are contained in the Supplement. A more extensive evaluation of the radiation, clouds, and aerosols in ECHAM6.3-HAM2.3 will be provided in the forthcoming papers of Tegen et al. (2018), Kokkola et al. (2018), Kühn et al. (2018), and Neubauer et al. (2018).

As a baseline for the following discussion, we use the ECHAM6.3-HAM2.3 simulation of Tegen et al. (2018) and an ECHAM6.1-HAM2.2 simulation from Zhang et al. (2014). The setups of these simulations differ from the ECHAM-HAMMOZ reference simulation described in this study. Not only were the ECHAM-HAM model versions run without the gas-phase module MOZ1.0, but a different setup for the nudging was also applied in both cases (temperature was not nudged for reasons described in Zhang et al. (2014); the nudging of surface pressure, vorticity, and divergence is similar to that of our ECHAM-HAMMOZ simulations). While the ECHAM-HAM simulations were specifically tuned to bring the TOA shortwave and longwave fluxes into the observed range (Hourdin et al., 2017), no special tuning was performed for the ECHAM-HAMMOZ base simulation. Furthermore, the ECHAM6.3-HAM2.3 simulation uses a different parameterization for the emission of sea salt aerosol particles (Long et al., 2011) than ECHAM6.1-HAM2.2 and the ECHAM-HAMMOZ reference simulation, which both use the parameterization from Guelle et al. (2001). Finally, the simulation period of ECHAM6.1-HAM2.2 covers the years 2000 to 2009, whereas ECHAM6.3-HAM2.3 and ECHAM-HAMMOZ were run from 2003 to 2012. Therefore, the purpose of the following comparisons is not specifically to evaluate the impact of gas-phase chemistry on aerosols, clouds, and radiation, but rather to provide a “sanity check” for ECHAM-HAMMOZ and place the simulations that are described in this study in context with the HAM literature.

Table 8 Global annual mean top of the atmosphere (TOA) energy budget, cloud-related properties, and aerosol mass burdens of the ECHAM-HAMMOZ base simulation, the ECHAM6.3-HAM2.1 (E63H23) and ECHAM6.1-HAM2.2 (E61H22) climate simulations of Tegen et al. (2018), and observations or multi-model mean values.

^a Taken from Fig. 1 of Wild et al. (2013). ^b Taken from Fig. 7.7 in Boucher et al. (2013); see references therein. ^c Johnson et al. (2016). ^d Stubenrauch et al. (2013). ^e Elsaesser et al. (2017). ^f Taken from Fig. 2 of Li et al. (2012). ^g Average of Table S1 of von der Haar et al. (2012). ^h Taken from Table 10 of Textor et al. (2006).

Global 10-year mean values of the top of the atmosphere (TOA) energy budget, cloud-related properties, and aerosol mass burdens are shown in Table 8 for the ECHAM-HAMMOZ reference simulation, the ECHAM6.1-HAM2.2 and ECHAM6.3-HAM2.3 simulations, and for a blend of observations and AeroCom multi-model mean values. The table footnotes provide references for each comparison. Maps showing the spatial distribution of aerosol burden can be found in Fig. S3.6 of the Supplement, and time series of aerosol burden are presented in Fig. S3.7.

The shortwave (SW) and longwave (LW) cloud radiative effects (CREs) are weaker in the ECHAM-HAMMOZ simulation than in ECHAM6.3-HAM2.3 and in the observations. The differences in SWCRE (LWCRE) are −6 W m⁻² (−2 W m⁻²) with respect to ECHAM6.3-HAM2.3 and −5 W m⁻² (−5 W m⁻²) with respect to the observations. These differences can be explained by a lower cloud cover and a smaller liquid water path in the ECHAM-HAMMOZ simulation. They are similar to differences observed between different ECHAM6.1-HAM2.2 simulations described in Zhang et al. (2014), in which different setups of the nudging were tested and compared to free-running simulations constrained only by sea surface temperature and sea ice cover. For ECHAM6.1-HAM2.2 they found a decrease in SWCRE (LWCRE) by about 5 W m⁻² (1 W m⁻²) together with a decrease in the liquid water path by about 10 g m⁻² when temperature was included in the nudging parameters. Additional impacts of nudging were seen in the water vapor path and in convective clouds. It is therefore likely that the differences in cloud and radiation fields between ECHAM-HAMMOZ and ECHAM-HAM can be explained by the different nudging setups used, although the different model tuning may also play a role here.

The differences in SWCRE and LWCRE lead to larger net shortwave TOA fluxes and an imbalance of the radiative TOA fluxes of 4.6 W m⁻², which is somewhat larger than the imbalances found in the two ECHAM-HAM simulations or in Johnson et al. (2016).

The 10-year global means of the aerosol mass burdens of the ECHAM-HAMMOZ simulations are comparable with the ECHAM-HAM simulations and AeroCom multi-model mean values of Textor et al. (2006). The particulate organic matter (POM) burdens in all ECHAM-HAM(MOZ) simulations are lower than in the AeroCom multi-model mean, which could be due to the simplistic treatment of secondary organic aerosols (Zhang et al., 2012). The POM burden in the base run is even lower than in the two comparison simulations. The reasons for this are not clear at present. Compared to ECHAM6.3-HAM2.3 the sea salt burdens in ECHAM-HAMMOZ and ECHAM6.1-HAM2.2 are substantially larger. This difference results from changing the sea salt emission scheme from Guelle et al. (2001) to Long et al. (2011) as described above. The multi-model mean from AeroCom (Textor et al., 2006) falls in between the two ECHAM-HAM versions.

Overall, the 10-year mean aerosol burden maps (Fig. S3.5) show a reasonable comparison between the patterns of ECHAM-HAMMOZ and ECHAM6.3HAM2.3. The gas-phase module MOZ1.0 in ECHAM-HAMMOZ does not interact directly with black carbon, particulate organic matter, sea salt, and mineral dust. However, the formation of sulfate aerosol is different (gas-phase H₂SO₄ is explicitly included in MOZ1.0, whereas it is parameterized in HAM) so that the largest differences between ECHAM-HAMMOZ and ECHAM6-HAM2.3 are expected for sulfate aerosol. Indeed, the sulfate burden (first row of Fig. S3.5) is lower everywhere in ECHAM-HAMMOZ than in ECHAM6.3-HAM2.3, but is quite comparable to ECHAM6.1-HAM2.2. The differences between the different ECHAM-HAM(MOZ) versions are in the range of the AeroCom models. The spatial distributions of black carbon, particulate organic matter, and mineral dust burdens are quite similar in ECHAM-HAMMOZ compared to ECHAM6.3-HAM2.3. For the sea salt burden, the spatial distribution of ECHAM-HAMMOZ agrees well with the one of ECHAM6.1-HAM2.2. The seasonal cycles of the global means of aerosol mass burdens in ECHAM-HAMMOZ are very similar to the ones in ECHAM6.3-HAM2.3 and ECHAM6.1-HAM2.2 (Fig. S3.6 in the Supplement). The differences are consistent with those discussed above.

Overall, it can be concluded that the nudged ECHAM-HAMMOZ base simulation produces reasonable results for TOA radiative fluxes, cloud-related properties, and aerosol mass burdens. All parameters are in the range of other nudged ECHAM-HAM simulations, AeroCom multi-model mean values, and observations. Differences between ECHAM-HAMMOZ and ECHAM6.3-HAM2.3 can be explained primarily by the different model setup, i.e., the use of a different sea salt emissions scheme (see also Tegen et al., 2018), and different settings for the dynamical nudging (Zhang et al., 2014). The largest direct impact of using the gas-phase module MOZ1.0 that could be identified is a 21 % decrease in sulfate aerosol mass burden, which is within the range of AeroCom models. A more thorough evaluation of the impact of using MOZ1.0 could be achieved by comparing model configurations with and without MOZ1.0 using an otherwise identical model setup, but this is beyond the scope of this study.

Shorter test simulations with the SALSA microphysics package reveal minor differences compared to the M7 simulations presented in this work (Stadtler et al., 2018).

7 Conclusions

ECHAM-HAMMOZ in its released version ECHAM6.3-HAM2.3-MOZ1.0 is a state-of-the-art chemistry–climate model with a comprehensive tropospheric and stratospheric chemistry package and two options to model aerosol processes with either a modal or a bin scheme.

A 10-year simulation from 2003 to 2012 in the default configuration was performed and has been evaluated with various observational data and compared to other model studies. The focus of the evaluation was placed on the year 2008. The model reproduces many of the observed features of total column ozone, polar stratospheric processes, tropospheric and surface ozone, and column and surface CO. Like many other models, ECHAM-HAMMOZ shows a high bias of surface ozone concentrations, but this bias is relatively modest. Global budgets of ozone and OH are in line with estimates from multi-model intercomparison studies and two individual models using either a similar chemistry scheme (CAM-CHEM) or a similar climate model (EMAC) as ECHAM-HAMMOZ.

The evaluation of the model run in the default configuration revealed two issues with respect to the global NO_x budget: (1) lightning emissions are only about 1.2 Tg(N) yr⁻¹ and thus a factor of 2 lower than the lower limit that is generally accepted by the community; and (2) the parameterizations of heterogeneous reactions constitute a too-strong sink of HNO₃. The aerosol model of ECHAM-HAMMOZ does not include an explicit treatment of nitrate, and therefore the reevaporation of HNO₃ that is lost to the aerosol phase does not occur. Three sensitivity simulations were performed, which corrected these issues. Unfortunately, the results from these simulations tend to increase model biases, and they particularly invigorate the tropospheric ozone chemistry and decrease the lifetime of CH₄. These changes occur almost exclusively in the tropics and may be related to issues with tropical dynamics in the ECHAM model, which also show up as precipitation bias. The evaluation of cloud and radiation budgets hints towards a possible radiative imbalance induced by the nudging setup. However, the precipitation bias has also been found in other simulations with ECHAM6, and the excessive ozone chemistry and OH concentrations are also a feature of EMAC, which is based on an earlier version of ECHAM (with modifications).