Earth System Chemistry integrated Modelling (ESCiMo) with the Modular Earth Submodel System (MESSy) version 2.51

Abstract

Three types of reference simulations, as recommended by the Chemistry–Climate Model Initiative (CCMI), have been performed with version 2.51 of the European Centre for Medium-Range Weather Forecasts – Hamburg (ECHAM)/Modular Earth Submodel System (MESSy) Atmospheric Chemistry (EMAC) model: hindcast simulations (1950–2011), hindcast simulations with specified dynamics (1979–2013), i.e. nudged towards ERA-Interim reanalysis data, and combined hindcast and projection simulations (1950–2100). The manuscript summarizes the updates of the model system and details the different model set-ups used, including the on-line calculated diagnostics. Simulations have been performed with two different nudging set-ups, with and without interactive tropospheric aerosol, and with and without a coupled ocean model. Two different vertical resolutions have been applied. The on-line calculated sources and sinks of reactive species are quantified and a first evaluation of the simulation results from a global perspective is provided as a quality check of the data. The focus is on the intercomparison of the different model set-ups. The simulation data will become publicly available via CCMI and the Climate and Environmental Retrieval and Archive (CERA) database of the German Climate Computing Centre (DKRZ). This manuscript is intended to serve as an extensive reference for further analyses of the Earth System Chemistry integrated Modelling (ESCiMo) simulations.

Introduction

The study of chemistry–climate interactions represents an important and, at the same time, difficult task of global change research. The emerging issues of climate change, ozone depletion and air quality, which are challenging from both, scientific and policy perspectives, are represented in chemistry–climate models (CCMs). Understanding how the chemistry and composition of the atmosphere may change over the 21st century is essential for preparing adequate adaptive responses or establishing mitigation strategies. The distribution and development of aerosols and reactive greenhouse gases is controlled by primary emissions, atmospheric chemistry, and physics including transport of air masses integrated over global scales. Projections of future climate change are coupled with changes in atmospheric composition, whose impacts extend from air quality to stratospheric ozone. Furthermore, chemically active species in the troposphere are more amenable to short-term manipulations by changes in emissions and are therefore of major policy relevance to both air quality and climate. Provision of high-quality, policy-relevant information on the current state of the climate and its possible future state, as well as options for adaptation, are strongly dependent on progress in this area.

Increasingly, the chemistry and dynamics of the stratosphere and troposphere are being studied and modelled as a single entity in global models. The European Centre for Medium-Range Weather Forecasts – Hamburg (ECHAM)/Modular Earth Submodel System (MESSy) Atmospheric Chemistry (EMAC) model was one of the first community models with this capability . For the first time, some of the Earth system models (ESMs) with interactive oceans participating in the Coupled Model Intercomparison Project Phase 5 CMIP5; had interactive chemistry . The World Meteorological Organization – United Nations Environment Programme (WMO/UNEP) “Scientific Assessment of Ozone Depletion: 2010” also featured several stratospheric models that included tropospheric chemistry, and one model with a coupled ocean. It was also a main recommendation of the Stratosphere–troposphere Processes And their Role in Climate (SPARC) Chemistry–Climate Model Validation Activity (CCMVal) report , that stratosphere-resolving CCMs should continue to evolve towards more comprehensive, self-consistent stratosphere–troposphere CCMs. These developments provide a pathway for including better representation of stratosphere–troposphere, chemistry–climate, and atmosphere–ocean coupling in CCMs and ESMs used for more robust predictions of future stratospheric ozone, climate change, and mutual influences .

Within the Earth System Chemistry integrated Modelling (ESCiMo) initiative chemistry–climate simulations have been conducted by the MESSy (http://www.messy-interface.org) Consortium with the EMAC model for special topics related to the national project of the DFG-Forschergruppe SHARP (Stratospheric Change and its Role for Climate Prediction) and the International Global Atmospheric Chemistry (IGAC)/SPARC Chemistry–Climate Model Initiative (CCMI; http://www.met.reading.ac.uk/ccmi/). These simulations have been carried out in support of upcoming WMO/UNEP ozone and IPCC climate assessments and will help to answer emerging science questions as well as to improve process understanding.

In this manuscript Sect.  documents briefly (mainly for the users of it) the updates of the MESSy and the EMAC model since , Sect.  describes the model set-ups, and Sect.  highlights some first analyses. Section  provides a summary and some first conclusions, followed by a section providing information about the code and data availability. Last but not least, Appendix  lists the applied on-line diagnostics. Extensive supplementary information is available in the Supplement.

New model developments

MESSy is a software package providing a framework for a standardized, bottom-up implementation of Earth system models with flexible complexity. “Bottom-up” means, the MESSy software provides an infrastructure with generalized interfaces for the standardized control and interconnection (coupling) of low-level ESM components (i.e. dynamic cores, physical parameterizations, chemistry packages, diagnostics), which are called submodels. MESSy comprises currently about 60 submodels (i.e. coded according to the MESSy standards) in different categories: infrastructure (i.e. framework) submodels, atmospheric chemistry-related submodels, physics-related submodels, and diagnostic submodels. The EMAC model uses the MESSy to link multi-institutional computer codes to the core atmospheric model, i.e. the 5th generation European Centre Hamburg general circulation model (ECHAM5; ). Since the publication of , MESSy – including the EMAC model – has undergone several updates.

Updates to the model infrastructure

The ECHAM5 “nudging” routines for simulations with specified dynamics have been updated to enable the usage of nudging data in netCDF (http://www.unidata.ucar.edu/software/netcdf/) input format. The corresponding preprocessing is now performed with cdo (climate data operators; https://code.zmaw.de/projects/cdo).

New infrastructure submodels for the unified import of external data (IMPORT), and grid definitions and transformations (GRID) have been implemented . The rediscretization via NCREGRID, as described by , is now part of those. Based on this advanced infrastructure, the submodels OFFLEM, ONLEM both, and DRYDEP for off-line, on-line calculated emissions, and dry deposition, respectively, have been revised (keeping their functionality) and renamed OFFEMIS, ONEMIS, and DDEP, respectively.

Updates of atmospheric chemistry-related submodels

The submodel MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere), used to simulate the chemical kinetics, has been updated and further revised. The rate coefficients were updated to the latest recommendations of the Jet Propulsion Laboratory JPL; and other recent publications (see Sect. ).

The submodel JVAL, used to calculate photolysis rate coefficients, has been updated to the most recent version see Sect. .

Updates of physics-related submodels

The ocean subsystem as described by , consisting of the submodels MPIOM based on the General (ocean) Circulation Model MPIOM version 1.3.0; , HD Hydrological Discharge; based on the work of, and A2O atmosphere to ocean coupling; are now included in the new MESSy version 2.51 (see Sect. ). The submodel A2O is responsible for the exchange of information between ocean and atmosphere (and vice versa), while HD simulates the riverine fresh water input into the ocean for balancing the water mass.

The radiation submodel RAD4ALL of development cycle 1 of MESSy , which was a first modularized version of the ECHAM5 radiation scheme, has been completely refined, further modularized, and split into the MESSy submodels RAD (with sub-submodel FUBRAD; see next item) for radiation calculations, AEROPT for the calculation of aerosol optical properties (see Sect. ), CLOUDOPT for the calculation of cloud optical properties, and ORBIT for Earth orbit calculations. The technical documentation with application examples will be published elsewhere ().

To increase the spectral resolution of the ultraviolet–visible (UV–Vis) region of the solar spectrum, the single UV–Vis band of the submodel RAD can be substituted by the sub-submodel RAD-FUBRAD at pressures lower than 70 $\mathrm{hPa}$ . The updates of FUBRAD were motivated by the need for a consistent flux profile over the complete vertical model domain, that is necessary when an interactive ocean is coupled to EMAC. The updates comprise

• an increase of the spectral resolution of the Chappuis band from one in the original version to 6 or 57 in the updated version, increasing the spectral resolution from 49 to 55 or 106 bands;

• a consistent flux profile by using non-scaled fluxes in the Chappuis band and also the usage of integrated fluxes in the non-absorbing band.

The surface processes of the ECHAM5 base model have been restructured as MESSy submodel SURFACE. The restructuring comprises the calculation of heat and water budgets of the surface (surf.f90)

Names in parentheses refer to the original ECHAM5 subroutine names.

, and lake temperatures and ice thicknesses (lake.f90). The ground temperature evolution and the temperature profile within the soil are estimated from the thermal diffusion equation (soiltemp.f90). Lake-ice temperature (licetemp.f90) and sea-ice temperature (sicetemp.f90) are calculated prognostically; prognostic calculation of the ice temperature (icetemp.f90) is optional and has been applied for all simulations presented here. The concepts and physics of the submodel SURFACE are described in detail by .

In summary, the physics-related submodels (i.e. the “E” in EMAC) have been mostly derived from the physics routines of the ECHAM5 base model , some have been further developed and include advanced or alternative parameterizations. Thus, ECHAM5 physics in EMAC is currently represented by the submodels AEROPT, CLOUDOPT, ORBIT, RAD, SURFACE (see above), GWAVE non-orographic gravity waves;, CONVECT convection;, and CLOUD. Convective tracer transport is simulated by the submodel CVTRANS . Vertical diffusion, orographic gravity waves, and large-scale advection flux-form semi-Lagrangian; are treated in EMAC by ECHAM5. The feedback of atmospheric chemistry to the hydrological cycle is controlled by the submodel H2O .

New diagnostic submodels

New diagnostic submodels have been added:

• SCALC see also Appendix  for an example usage; is used for Simple CALCulations with “channel objects” via namelist. In the current MESSy release 2.51 summation and scaling of channel objects is implemented. The term channel object was introduced as part of the MESSy terminology by . In brief, it describes a specific Fortran95 structure comprising the data and corresponding meta-data of prognostic and diagnostic variables according to an object-oriented approach. The individual model components (i.e. what we call submodels) operate on these channel objects. SCALC, in particular, is used to provide, defined by namelist, new channel objects (e.g. the total loss rate of a reactive compound), consisting of the sum of (optionally scaled) individual objects (e.g. the process specific loss rates of that compound).

• The submodel TBUDGET (see Appendix ) is used to calculate budgets of reactive compounds.

• The submodel CONTRAIL (see Appendix ) is used to calculate the potential contrail coverage and the potential contrail cirrus coverage.

Overview of ESCiMo reference simulations. The CCMI notation Ref-C1/2 is abbreviated with RC1/2, SD denotes specified dynamics. Different model configurations are indicated by -base-, -aero-, -aecl- and -oce-. The two digit numbers indicate the specific simulations (e.g. to distinguish the vertical resolution and specific SD set-up). An appended letter “a” indicates a sensitivity study (see remarks and Sect. 3.12). For detailed explanations, see text. The line colours refer to the figures below.

[Figure omitted. See PDF]

Simulation configurations and set-ups

General aspects of the model set-ups

Three different sets of reference simulations have been suggested by CCMI , namely free-running hindcast simulations from 1960 to 2010 (REF-C1), hindcast simulations with specified dynamics from 1980 to 2010 (REF-C1SD), and combined free-running hindcast and projection simulations from 1960 to 2100 (REF-C2). With EMAC we conducted these simulations in a T42 (triangular) spectral resolution of the ECHAM5 base model (corresponding to a quadratic Gaussian grid of 2.8${}^{\circ }$ $\times$ 2.8${}^{\circ }$ in latitude and longitude) in different flavours as listed in Tables  and . We simulated with two different vertical resolutions, once with 90 (L90MA) model levels and once with 47 (L47MA) model levels, both reaching up to 0.01 $\mathrm{hPa}$ (mid- of uppermost layer) into the middle atmosphere (MA; approximately 80 $\mathrm{km}$). The vertical layer structures of the hybrid pressure levels ($p(i,j,k,t)=h_a\left(k\right)+h_b\left(k\right)\cdot p_{\mathrm{s}}(i,j,t)$, where $p$ is pressure, $i,j$ the horizonal and $k$ the vertical grid indices, $t$ time, $p_{\mathrm{s}}$ the surface pressure, and $h_a$ and $h_b$ the hybrid coefficients) are visualized in Fig. S4

Figures and tables named S$n$ (with $n$ $=$ 1, 2, 3, …) referred to in the text are in the document ESCiMo_supplement.pdf. Those named E$n$ (with $n$ $=$ 1, 2, 3, …) are in the document ESCiMo_emissions.pdf. Both documents are part of the Supplement.

. The time step length was 720 $\mathrm{s}$ in the T42L90MA and 600 $\mathrm{s}$ in the T42L47MA simulations, respectively. For the T42L47MA resolution, we had to optimize several model parameters as described in further detail in Sect. .

All simulations (except for those with specified dynamics, SD) start in January 1950, to have a 10-year-long spin-up period (1950–1959), initialized from already spun-up states of previous simulations; see Sect. . The simulation with coupled interactive ocean (RC2-oce-01) was spun-up in a two-stage procedure over 500 years in total (see Sect.  for details).

The hindcast simulations with specified dynamics (SD) have been branched off from restart files (1 January 1979) of the corresponding free-running hindcast simulations and “nudged” by Newtonian relaxation towards ERA-Interim reanalysis data , which are available with a 6-hourly time resolution from the year 1979 onwards. The Newtonian relaxation (nudging) of the ECHAM5 base model is applied in spectral space for the prognostic variables divergence, vorticity, temperature, and the (logarithm of the) surface pressure. The corresponding relaxation times applied were 48, 6, 24, and 24 $\mathrm{h}$, respectively. However, the nudging strengths are not applied homogeneously in the vertical: the boundary layer and the stratosphere–middle atmosphere above 10 $\mathrm{hPa}$ are not nudged with transition layers of intermediate strengths in between. The vertical profiles of the relative nudging strengths for both vertical resolutions are displayed in Fig. S5. The nudging further implies that the sea surface temperatures (SSTs) and the sea-ice concentrations (SICs) are used from ERA-Interim reanalysis data, whereas the free-running simulations (RC1 and RC2) are forced by other external SSTs and SICs (see Sect. ). In addition to the RC2-base simulations (free-running hindcast and projection), which are forced by prescribed SSTs/SICs (see Sect. ), we conducted a simulation with an interactively coupled ocean model . The details are described in Sect. .

As listed in Table , the SD simulations have been performed twice (in both vertical resolutions) with two different settings: in two cases the “wave zero” (i.e. the global mean) temperature ($\stackrel{\mathrm{‾}}{T}$) was included for the Newtonian relaxation (RC1SD-base-07 in T42L90MA and RC1SD-base-08 in T42L47MA), in two other cases (RC1SD-base-10 in T42L90MA and RC1SD-base-09 in T42L47MA) it was omitted. The idea was to correct (or not) potential temperature biases and to investigate the response of the model, in particular the chemical state of the model atmosphere.

For investigations about the role of tropospheric aerosol for chemistry, we additionally performed free-running hindcast simulations with interactively calculated aerosol replacing the prescribed aerosol of the RC1-base simulations, once without (RC1-aero-06/07) and once with effect on the clouds (RC1-aecl-01/02). Further details are given in Sect.  and . Last, but not least, we performed additional sensitivity simulations (those appended by letter “a” in Table ) for reasons outlined in Sect. .

Data from the model were mostly output as 10-hourly global snapshots or averages. All simulation set-ups were equipped with extended on-line diagnostics as detailed in Appendix . The high complexity of the applied model set-ups and the total number of simulation years required a considerable amount of computational resources and led to almost 2 $\mathrm{PetaByte}$ of primary model output (Table ).

Parameter optimization for T42L47MA

In preparation of the simulations with interactive chemistry a number of 10- to 20-year test simulations have been performed with EMAC at a resolution of T42L47MA to determine optimal parameter settings. In these simulations the interactive chemistry has been switched off, whereas the same set-up as for the chemistry simulations has been used for convection submodel CONVECT;, cloud cover submodel CLOUD;, and non-orographic gravity waves submodel GWAVE;. Optimization of parameters is necessary for the convection parameterization to get a balanced radiation budget at the top of the atmosphere (TOA), for the quasi-biennial oscillation (QBO) nudging to get a realistic QBO amplitude for the L47MA resolution, and for the gravity wave parameterization, where the choice of the parameter rmscon (the root mean square of the gravity-wave-induced horizontal wind speed at the launch level in $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) has an impact on the polar vortex strength.

The test simulations are performed to achieve a global, annual average equality of the net incoming shortwave (SW) radiation with the outgoing longwave (LW) radiation at the uppermost model level (i.e. TOA). The test simulation was performed under conditions for the year 2000 for green house gases (GHGs), ozone depleting substances (ODSs), SSTs and SICs (10-year average of the Hadley Centre Sea Ice and Sea Surface Temperature data set (HADISST) monthly SSTs and SICs between 1995 and 2004). Comparing the TOA balance of the L47 simulations with interactive chemistry reveals an annual, global average from 1995 to 2004 of $-$0.26 and 0.41 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for RC1-base-08 and RC2-base-05, respectively.

Comparison to the test simulation without interactive chemistry (0.1 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$) shows that these values are still in the range of $\pm$0.5 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, only slightly larger than the uncertainty range from observations. give an estimate for the TOA radiation balance of 0.6 ($\pm$0.4) $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ for the decade 2000–2010 derived from satellite observations.

Used computational resources and size of simulated data

Simulation	Resolution	Simulated years		Wall clock	Data size${}^{\text{a}}$	CPU h${}^{\text{b}}$	Data size${}^{\text{a}}$
				[days year${}^{-\mathrm{1}}$]	[TByte year${}^{-\mathrm{1}}$]		[TByte]
RC1-base-07	T42L90MA	1950–2011	(62)	2.58	2.54	491 444	157.79
RC1-base-07a	T42L90MA	1990–2010	(21)	2.58	2.55	166 597	53.52
RC1-base-08	T42L47MA	1950–2011	(62)	1.81	1.35	344 194	83.68
RC1-base-08a	T42L47MA	1990–2010	(21)	1.82	1.35	117 448	28.37
RC1-aero-06${}^{\text{c}}$	T42L90MA	1950–1998	(48)	3.84	3.74	573 680	182.62
RC1-aero-07${}^{\text{c}}$	T42L90MA	1990–2011	(22)	3.84	3.74	259 666	82.56
RC1-aecl-01${}^{\text{d}}$	T42L90MA	1958–1972	(14)	3.88	3.87	175 125	57.63
RC1-aecl-02${}^{\text{d}}$	T42L90MA	1965–2011	(47)	3.88	3.87	561 085	182.62
RC1SD-base-07	T42L90MA	1979–2013	(35)	2.74	2.60	295 070	90.86
RC1SD-base-08	T42L47MA	1979–2013	(35)	1.96	1.34	210 639	48.24
RC1SD-base-09	T42L47MA	1979–2013	(35)	1.96	1.38	210 792	48.26
RC1SD-base-10	T42L90MA	1979–2013	(35)	2.70	2.60	290 347	90.97
RC1SD-base-10a	T42L90MA	2000–2014	(15)	2.76	2.61	127 347	39.11
RC2-base-04	T42L90MA	1950–2099	(150)	2.57	2.55	1 183 680	382.23
RC2-base-05	T42L47MA	1950–2099	(150)	1.81	1.35	833 640	202.72
RC2-oce-01	T42L47MA/GR30L40	1950–2100	(151)	1.86	1.35	861 610	204.40
SUM						6 702 364	1935.58

${}^{\text{a}}$ Including restart files (every 3 months and QTIMER triggered at the end of the scheduler wall-clock limit); ${}^{\text{b}}$ on 4 nodes of an IBM Power6 in SMTP mode (i.e. with 64 tasks/node); ${}^{\text{c}}$ temporal overlap of 06 with 07: January 1990–August 1998; ${}^{\text{d}}$ temporal overlap of 01 with 02: January 1965–October 1972.

Balancing the radiation budget

To achieve a balanced radiation budget at the TOA, the tuning of cloud parameters is commonly applied e.g.. The default settings for the T42L47MA set-up of EMAC result in a TOA radiation imbalance of $-$3.3 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, i.e.the outgoing LW radiation (OLR) is larger than the net SW radiation. The prescribed SSTs limit the ability of the model to adapt to a radiation imbalance with increasing/decreasing SSTs, which would also influence convection and clouds, and thereby possibly balance the radiation budget by increasing/decreasing the outgoing LW radiation. With this kind of model set-up the negative imbalance can be reduced by increasing the net SW radiation through changes of the parameters that affect cloud properties, thus reducing the planetary albedo. There are a number of parameters in the Tiedtke convection parameterization, which have an influence on the net SW radiation, as summarised by . In a series of test simulations optimized values of the following parameters have been identified that lead to an almost balanced radiation budget at the TOA: the relative convective cloud mass flux above level of non-buoyancy (cmfctop $=$ 0.35, default: 0.30, possible range: 0.10–0.38) and the entrainment rate for deep convection (entrpen $=$ $0.5\times {10}^{-\mathrm{4}}$ ${\mathrm{m}}^{-\mathrm{1}}$, default: $1.0\times {10}^{-\mathrm{4}}$ ${\mathrm{m}}^{-\mathrm{1}}$, possible range: 0.3–$3.0\times {10}^{-\mathrm{4}}$ ${\mathrm{m}}^{-\mathrm{1}}$). The values for the possible ranges of these two parameters have been tested by . Both changes, increasing the parameter cmfctop and decreasing the parameter entrpen, increase the net SW radiation at TOA, as the low-level clouds tend to be less and thinner. Applying both altered parameters in combination leads to an 4.3 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ increase of the net SW radiation, and an 0.9 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ increase of the OLR at TOA, which results in an almost balanced TOA radiation budget of 0.1 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$. For the T42L90MA simulations, the default values have been used (see Table S2). For the simulation with coupled ocean model, cloud parameters have also been modified (see Sect.  and Table S2).

QBO nudging

The vertical resolution of L47MA is not sufficient to generate the quasi-biennial oscillation (QBO) of the zonal winds in the lower equatorial stratosphere internally. Therefore, the zonal winds near the Equator are relaxed (i.e. nudged with submodel QBO) towards a zonal mean field with a Gaussian profile in the latitudinal direction, which has been derived from the observed zonal mean zonal winds near the Equator, to get the correct amplitude and phase of the observed QBO (see Fig. S7). The nudging is applied in the altitude range between 10–90 $\mathrm{hPa}$, with full nudging weights (i.e. 1.0) from 20 to 50 $\mathrm{hPa}$, levelling off to 0.3 (0.2) at the upper (lower) edge of the nudging region . The latitudinal range is confined to 12.6${}^{\circ }$ S–12.6${}^{\circ }$ N, with full nudging at latitudes from 7${}^{\circ }$ S–7${}^{\circ }$ N. The relaxation timescale is set to 10 $\mathrm{days}$. Although the internal generation of a QBO is a feature of the L90MA set-up of EMAC , the zonal winds near the Equator are slightly nudged with a relaxation timescale of 58 $\mathrm{days}$ in all EMAC simulations at T42L90MA resolution, to get the correct phasing of the observed QBO.

Differences of the daily mean climatology (19 $\mathrm{years}$) between a T42L47MA EMAC simulation with the gravity wave parameter rmscon set to 0.92 and a simulation where it is set to 0.96; for the zonal mean zonal winds at 60${}^{\circ }$ S in $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ (top); for the zonal mean temperatures averaged from 71.2 to 87.9${}^{\circ }$ S in $\mathrm{K}$ (bottom). The shading indicates the differences that are significant at the 95 $\mathrm{\%}$ level, estimated with a Student's $t$ test.

[Figure omitted. See PDF]

Tuning of gravity waves

The strength of the momentum deposition in the stratosphere and mesosphere by non-orographic gravity waves, thought to be released at the launch level, which is chosen to be near a pressure of 643 $\mathrm{hPa}$, is controlled in submodel GWAVE by the parameter rmscon. The parameter is used in the calculation of the cutoff wave number in each model level, controlling the wave number at which the wave breakdown begins . Decreasing rmscon leads to a larger cutoff wave number and wave breaking occurs less often, which has the effect of reducing the disturbance of the polar vortex by deposited momentum flux in the stratosphere and mesosphere. In earlier EMAC simulations at T42L90MA resolution this parameter was chosen to be 0.96, which often has led to a relatively weak Antarctic polar vortex with a warm bias in the stratosphere. The setting of rmscon to 0.92, for both T42L47MA and T42L90MA resolutions, was found to be optimal for the strength of the Antarctic polar vortex and was determined by a series of test simulations with rmscon set to 0.90, 0.92, 0.96, and 0.98. As shown for the daily mean climatology of the zonal wind at 60${}^{\circ }$ S and the temperatures for a polar cap average from 71.2–87.9${}^{\circ }$ S (Fig. ), this choice of the parameter leads to a strengthened Antarctic polar vortex and a colder stratosphere, especially during August and September.

The effects on the polar vortex strength in the Northern Hemisphere are more diverse, with alternating periods of an intensified and attenuated polar vortex. This is also reflected in the temperature changes of the polar cap average from 71.2 to 87.9${}^{\circ }$ N. Most of these differences are, however, not statistically significant and therefore not shown.

Sea surface temperatures and sea-ice concentrations

SSTs and SICs for the RC1 simulations are prescribed following the global data set HadISST provided by the UK Met Office Hadley Centre available via http://www.metoffice.gov.uk/hadobs/hadisst/;. The data set is based on merged satellite and in situ observations. For the RC1SD simulations, SSTs and SICs as given by ERA-Interim were used for consistency with the nudging (see Sect. ). In the global mean, the HadISST and ERA-Interim SSTs and SICs are almost identical (Fig. ). The distribution of SSTs is almost identical in ERA-Interim and HadISST as well (Figs. S8 and S9), and SICs differ slightly in the pattern by up to $\pm$20 $\mathrm{\%}$ (Fig. S12).

SSTs and SICs for simulations extending into the future (i.e. RC2) are taken from simulations with the Hadley Centre Global Environment Model version 2 - Earth System (HadGEM2-ES) Model . These simulations were performed for the CMIP5 multi-model data sets, and have been made available via the CMIP5 data archive at PCMDI (Program for Climate Model Diagnosis and Intercomparison; available at http://www-pcmdi.llnl.gov). For years up to 2005, the “historical” simulation with HadGEM2-ES is used. Afterwards, the “RCP6.0” simulation, which is initialized with the historical simulation, is employed. Details of those simulations are described by . The simulation of SSTs in HadGEM2 has been significantly improved over the predecessor model HadGEM1, namely, the prominent cold bias in HadGEM1 has been reduced in HadGEM2, and the representation of ENSO has been improved see also Figs.  and S10;. The sea-ice extent has been well reproduced, remaining within 20 $\mathrm{\%}$ of observed values during most of the year see also Figs.  and S13;, even though we find deviations of up to 80 $\mathrm{\%}$ locally.

It would have been desirable to take future SSTs and SICs from simulations with a climate model based on the same atmospheric base model as EMAC. However, the suitable simulations following RCP6.0 were not performed with the corresponding model (Max Planck Institute Earth system model, MPI-ESM).

Transient spectrally resolved irradiances

To account for the solar variability, daily spectrally resolved irradiances (SSI) from the Naval Research Laboratory Solar Spectral Irradiance Model (NRLSSI) and the daily total solar irradiance TSI; have been used. Data have been prepared as input for the radiation (RAD) sub-submodel FUBRAD , here applied with 55 spectral bands, and for the photolysis calculations in the submodel JVAL (see Sect. ). Corresponding time series are shown in Figs. S14 and S15. Note that the temporal evolution of the photon fluxes derived for JVAL (not shown) directly follow the shown spectral irradiances.

The future solar forcing, to be used for the projections, has been prepared according to the solar forcing used for CMIP5 simulation of HadGEM2-ES, where the SSTs and SICs are taken from  (2011; see also Sect. ). It consists of repetitions of an idealized solar cycle connected to the observed time series in July 2008. This has been applied consistently for all projections with prescribed SSTs (RC2-base) and for the simulations with interactively coupled ocean model (RC2-oce-01). Here, we deviate from the CCMI recommendations consisting of a sequence of the last four solar cycles (20–23).

Atmospheric chemistry set-ups

Gas-phase chemistry and photolysis

For the chemical kinetics, we have used the submodel MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere). Based on the code described by and , a revised version was used for the simulations. All rate coefficients were updated to the latest recommendations by JPL and other recent publications. Apart from these updated rate coefficients, the product distributions were also updated for a few reactions (e.g. ${\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{4}}$ + ${\mathrm{O}}_{\mathrm{3}}$). In addition, previously neglected, chemically inert or ubiquitous products like ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ and ${\mathrm{O}}_{\mathrm{2}}$, have now been added in order to fix the mass balance of some reactions. The most recent version of the submodel JVAL was used to calculate photolysis rate coefficients ($J$ values), as described by . Spectrally resolved time series of photon fluxes have been used consistently with the transient spectrally resolved irradiances used for the radiation calculation (see Sect. ).

For the base simulations (RC1-base, RC1SD-base, and RC2-base; see Sect.  and Table ), the chemical mechanism was selected with the batch file CCMI-base-02.bat (see Figs. S1 and S2). Briefly, the mechanism considers the basic gas-phase chemistry of ozone, methane, and odd nitrogen. Alkanes and alkenes are included up to ${\mathrm{C}}_{\mathrm{4}}$. Alkynes and aromatics are not considered in our mechanism. Halogen chemistry includes bromine and chlorine species. For the chemistry of isoprene plus a few selected non-methane hydrocarbons (NMHCs), we used version 1 of the Mainz Isoprene Mechanism (MIM1) based on . Heterogeneous reactions of dinitrogen pentoxide (${\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$), halogen nitrates (${\mathrm{ClNO}}_{\mathrm{3}}$, ${\mathrm{BrNO}}_{\mathrm{3}}$), and hypohalous acids ($\mathrm{HOCl}$, $\mathrm{HOBr}$) are also included. Since $\mathrm{Hg}$ chemistry is not considered in this study, all $\mathrm{Hg}$ reactions were switched off. In total, the mechanism is described by 310 reactions of 155 species. For the simulations with interactive tropospheric aerosol (RC1-aero and RC1-aecl; see Sect.  and Table ) the mechanism selected with the batch file CCMI-aero-02.bat (see Fig. S3) contains additional sulfur reactions (5 additional species and 11 additional reactions). A complete list of chemical reactions, rate coefficients, and references can be found in the file ESCiMo_MECCA_mechanism.pdf in the Supplement.

Aqueous-phase chemistry and wet deposition

Aqueous-phase chemistry in clouds and wet deposition are simulated with the help of the combined explicit scavenging submodel SCAV , which calculates the uptake/release to/from the gas and aqueous phase and subsequent wet deposition. In contrast to more simplified schemes, dissociation and aqueous-phase redox reactions are also explicitly calculated, e.g. the sulfur(IV) to sulfur(VI) oxidation, such that the effective exchange between gas and liquid phase is taken into account. The scheme also includes nitric acid (${\mathrm{HNO}}_{\mathrm{3}}$) uptake on ice particles (except for polar stratospheric cloud (PSC) particles) according to a Langmuir uptake and subsequent denitrification by sedimenting ice particles. Wet deposition is calculated from the in-cloud (and subsequent conversion of in-cloud to in-precipitation) and in-precipitation chemical concentrations for both large-scale and convective clouds. The chemical species and reactions, which comprise the liquid-phase chemical mechanism, can be found in the document ESCiMo_SCAV_mechanism.pdf, which is part of the Supplement.

Stratospheric heterogeneous chemistry and PSCs

The submodel MSBM (Multi-phase Stratospheric Box Model) simulates the number densities, mean radii, and surface areas of the sulfuric acid aerosols and the different polar stratospheric cloud (PSC) particles (supersaturated ternary solution (STS), nitric acid trihydrate (NAT), and ice). Further, the rate coefficients of all heterogeneous reactions (in the stratosphere) are calculated and used by the submodel MECCA (merged with the corresponding tropospheric values by its sub-submodel MECCA_KHET; see also Sect. ).

For the formation of NAT, a kinetic growth NAT parameterization is used with the assumption of a necessary supercooling of 3 $\mathrm{K}$ . Hence, homogeneous NAT nucleation is also possible. The formation of STS is based on , the formation of ice particles is based on the thermodynamic approach of . For the ice nucleation 20 $\mathrm{\%}$ ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ supersaturation is assumed to be necessary. A trapezoid scheme is used for the sedimentation of PSC particles.

Heterogeneous reaction rates and their temperature dependencies on NAT are calculated according to the parameterization of based on the measurements of . The heterogeneous reaction rate coefficients on liquid particles are taken from and . The uptake coefficients and reaction probabilities for ice particles are taken from .

Here, stratospheric ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ mixing ratios have been prescribed with a time series provided by the CCMI database (B. Luo, personal communication, 2013; ftp://iacftp.ethz.ch/pub_read/luo/ccmi/) for the period 1960 to 2011, with some artificial spikes removed in the years 1979–1982 (C. Brühl, personal communication, 2014). The corresponding values have been determined from data of several satellite instruments (SAGE, SAGE II, CALIPSO, GOMOS). The resulting time series is available as monthly zonal averages, 5${}^{\circ }$ in latitude and on 70 pressure levels between 530 and 3 $\mathrm{hPa}$. For the spin-up period (1950 to 1959), data from the year 1960 were used; for the years after 2011, data from the year 2011 were used.

Dry deposition and sedimentation

Dry deposition is an important sink for gas and aerosol-phase species. Additionally, sedimentation leads to a significant loss of aerosol particles from the atmosphere. Within the MESSy submodel DDEP formerly named DRYDEP; see Sect. ;, dry deposition velocities are calculated following the big leaf approach as proposed by .

The sedimentation of aerosol particles depends among others on the aerosol density and size. Aerosol particles are sedimented using the simple upwind scheme of the MESSy submodel SEDI .

In the simulations without prognostic aerosol chemical and microphysical properties (i.e. all except for -aero- and -aecl-), sedimentation fluxes are calculated by SEDI for the residual aerosols originating from evaporation of clouds and precipitation leading to particles. In these cases, particle size distribution (mean radius $=$ $\mathrm{5}\times {10}^{-07}$ $\mathrm{m}$, $\mathit{\sigma }=2.0$) and particle density ($\mathit{\rho }=1841.0$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$) are prescribed.

Initial conditions of trace gases

Initial conditions for January 1950 (start of the simulation spin-up phase) have been generated by scaling simulated atmospheric mixing ratios from end of December 2000 of a previous EMAC simulation . Short-lived reactive species, e.g. nitrogen oxides, nitric acid (${\mathrm{HNO}}_{\mathrm{3}}$), and ozone, were initialized directly with the mixing ratios from the year 2000. During the spin-up phase (1950–1959) these species undergo processing and the mixing ratios adjust on timescales of the order of several months. Species with longer atmospheric lifetime were initialized with scaled mixing ratios, as the adjustment during the simulation would require several years. Specifically, not only the greenhouse gases carbon dioxide (${\mathrm{CO}}_{\mathrm{2}}$), nitrous oxide (${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$), and methane (${\mathrm{CH}}_{\mathrm{4}}$), but also more reactive species, such as carbon monoxide ($\mathrm{CO}$) and the chlorofluorocarbons (CFCs), were initialized with mixing ratio distributions equivalent to reduced mixing ratios scaled from year 2000 conditions. Initial tracer fields have been generated through scaling with constant factors (see Table S1), which have been derived from the temporal evolution of the corresponding primary emissions in the period 1950 to 2000.

Prescribed boundary conditions: gas-phase species

For long-lived species, which are relevant to atmospheric chemistry and climate, pseudo-emissions are calculated by the submodel TNUDGE . This approach is chosen, due to most models' inability to correctly simulate the corresponding trends, if direct emissions are prescribed. This issue is in part related not only to uncertainties in emission estimates themselves, but also in the difficulties to accurately simulate the species' lifetimes. Therefore, with TNUDGE the simulated mixing ratios in the lowest model layer are relaxed by Newtonian relaxation to observed or projected surface mixing ratios. Species that are prescribed with TNUDGE are the greenhouse gases (${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ and ${\mathrm{CH}}_{\mathrm{4}}$), ozone depleting substances (CFCs: ${\mathrm{CFCl}}_{\mathrm{3}}$, ${\mathrm{CF}}_{\mathrm{2}}{\mathrm{Cl}}_{\mathrm{2}}$, ${\mathrm{CH}}_{\mathrm{3}}{\mathrm{CCl}}_{\mathrm{3}}$, ${\mathrm{CCl}}_{\mathrm{4}}$; HCFCs: ${\mathrm{CH}}_{\mathrm{3}}\mathrm{Cl}$, ${\mathrm{CH}}_{\mathrm{3}}\mathrm{Br}$; Halons: ${\mathrm{CF}}_{\mathrm{2}}\mathrm{ClBr}$, ${\mathrm{CF}}_{\mathrm{3}}\mathrm{Br}$), ${\mathrm{H}}_{\mathrm{2}}$, and ${\mathrm{SF}}_{\mathrm{6}}$ (Figs. E1–E6). In the RC1-aero and RC1-aecl simulations $\mathrm{COS}$ has been prescribed as well see Fig. E7.

For all species except $\mathrm{COS}$, surface mixing ratios from several globally distributed observation sites were taken from the Advanced Global Atmospheric Gases Experiment (AGAGE; http://agage.eas.gatech.edu) and the National Oceanic and Atmospheric Administration/Earth System Research Laboratory (NOAA/ESRL, http://www.esrl.noaa.gov). Where no observations were available, historical and projected mixing ratios of the greenhouse gases and ${\mathrm{SF}}_{\mathrm{6}}$ were taken from and extended with the RCP6.0 scenario as proposed by . Mixing ratios of ozone depleting substances were based on the halogen scenario A1 from Table 5-A3 of the .

As the mixing ratios recommended by CCMI are only available as global and annual averages, we conducted the following approach: a climatological annual cycle and the latitudinal distribution of each species are calculated from the observed mixing ratios (see above) and applied to the prescribed globally and annually averaged mixing ratios. Thus, we compiled for each species a monthly and latitudinally varying time series from 1950 to 2100. For ${\mathrm{H}}_{\mathrm{2}}$ we extrapolated observed mixing ratios linearly with a typical seasonal cycle to get a time series from 1950 to 2100.

Initially, we planned to prescribe the recommended time series only during the periods when no observed mixing ratios were available. However, the comparison of observed and recommended values (from CCMI) show partly significant differences (compare Figs. E1–E6), so we decided to use for the RC2 simulations only the recommended mixing ratios in order to apply consistent time series.

In the RC1 simulations, prescribed boundary conditions of the last available year (typically 2010 or 2011) have been used for later years as well, with one exception: in RC1SD-base-10a (see Table  and Sect. ) we used the RC2 set-up for the years 2012 and later.

Prescribed aerosol

In the standard simulations (-base-), the properties of the atmospheric aerosol (including volcanic aerosol) have been prescribed to take the interactions with radiation and heterogeneous chemistry into account. Only in the simulations RC1-aero and RC1-aecl (see Sect. ) prognostic aerosol is calculated with a modal scheme with four log-normal modes separated into hydrophilic internally mixed and hydrophobic externally mixed particles. Furthermore, in RC1-aecl the direct feedback of the aerosol on both chemistry and the circulation is explicitly taken into account. The prescribed aerosol effects are separated into the aerosol surface area, representing chemical effects via heterogeneous chemistry, and the radiative properties influencing the radiation budget.

Radiative properties

In the -base- simulations as well as in RC1-aero the radiative properties of the aerosol, which are used in the radiation calculation scheme of the model to estimate scattering and absorption by aerosol particles, have been calculated based on climatologies, i.e. on off-line prescribed data. The AEROPT submodel (for the calculation of aerosol optical properties) has been extended to allow for, in addition to the individual calculations, also the merging of two data sets for the optical aerosol properties. In the model set-ups applied here, the values for the tropospheric aerosol were calculated on-line and the values for the stratosphere have been determined from satellite data (see below). These two data sets are merged using tropospheric values below 500 $\mathrm{hPa}$, stratospheric values above 300 $\mathrm{hPa}$, and linearly interpolated values between 500 and 300 $\mathrm{hPa}$.

For the troposphere the Tanre climatology , as used in the standard calculations for EMAC (and one option of the extended AEROPT submodel), is applied to determine the radiative effects of aerosol particles. This low-resolution spectral aerosol climatology data use the actual relative humidity to determine the total aerosol extinction, single scattering albedo, and asymmetry factor, which are fed into the radiation calculation.

For the stratosphere, the aerosol radiative properties, i.e. the extinction, the single scattering albedo, and the asymmetry factor, were provided by the CCMI database, consistently derived from satellite observations as the stratospheric ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$ mixing ratio (monthly zonal averages, 5${}^{\circ }$ latitude resolution; see Sect. ). The corresponding values were inter-/extrapolated to the radiation bands of the EMAC model and provided with a 500 $\mathrm{m}$ grid spacing in the vertical. As for ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$, data from the year 1960 have been used for the spin-up period 1950 to 1959, and data from the year 2011 have been used for years 2011–2100.

Surface area for heterogeneous chemistry

The aerosol surface area is required for the calculation of heterogeneous chemical reactions on atmospheric particles. The most important reactions are the conversions of ${\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}$ into ${\mathrm{HNO}}_{\mathrm{3}}$ or aerosol nitrate. The required field for the aerosol surface area is a merged data set of values for the troposphere and the stratosphere.

The tropospheric aerosol surface-area concentration climatology (monthly averages repeated every year) for the base simulations has been derived from the LOW_AIR simulation by . This simulation was performed with EMAC (MESSy version 1.4) coupled to the aerosol submodel MADE. It covers a period of 10 $\mathrm{years}$ (1996–2005) in nudged mode, using the T42L19 resolution and a simplified chemical mechanism, including basic tropospheric reactions and the sulfur cycle. Emissions were based on the CMIP5 inventory for the year 2000 .

Based on these prescribed aerosol surface concentrations, the sub-submodel MECCA_KHET calculates the reaction coefficients for heterogeneous reactions in the troposphere and merges the latter with the corresponding stratospheric values as calculated by the submodel MSBM (see Sect. ).

Simulations with tropospheric aerosol and implications for the radiation budget

In contrast to the base scenarios interactive prognostic aerosol calculations are performed in the RC1-aero and RC1-aecl simulations. To determine the physical and chemical properties of the atmospheric aerosol, the EMAC aerosol submodel GMXE has been applied. GMXE calculates the microphysical properties of aerosols based on nucleation of sulfuric acid, condensation, and phase partitioning of semi-volatile inorganic components according to ISORROPIA-2 , coagulation, and a dynamical shift between the size categories on 4 log-normal modes. Additionally, aerosol particles are distinguished by their solubility, resulting in seven size categories (three externally mixed hydrophobic, and four internally mixed hydrophilic modes). The condensation routines and primary particle emissions also determine the overall chemical composition of the aerosol particles. Aerosol emissions as provided by the CCMI database are mapped on-line to the model grid using the algorithm of in the new submodel IMPORT , vertically redistributed via OFFEMIS , and assigned to the corresponding aerosol species within GMXE. The physical loss processes sedimentation, dry deposition, wet scavenging, and aerosol activation are considered in the corresponding EMAC submodels .

The RC1-aero and RC1-aecl simulations consider the effect of phase partitioning of ${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}$, ${\mathrm{HNO}}_{\mathrm{3}}$, $\mathrm{HCl}$, and ${\mathrm{NH}}_{\mathrm{3}}$, as well as the interaction of these compounds with primary aerosol species such as ${\mathrm{Na}}^{+}$, organic carbon (OC), black carbon (BC), and dust, resulting in overall feedbacks on the chemical composition of the atmosphere. Additionally, the interactively simulated aerosol also provides a consistent aerosol surface for the troposphere, which is used for the above-mentioned heterogeneous chemistry on aerosol particles.

Note, that there is no feedback between the aerosol properties and the radiation in the RC1-aero scenario, such that a configuration identical to the base scenarios is applied. Nevertheless, the optical properties of the interactive aerosol are determined diagnostically, i.e. not used in the radiation calculation. Instead, this simulation focuses on the interactions of the atmospheric aerosol and the gas-phase chemistry.

The RC1-aecl simulation, on the other hand, uses the interactively calculated aerosol also to determine the optical properties with the help of the AEROPT submodel , replacing the climatology in the troposphere. Note that the extensions of have not been applied here and that for the stratosphere, the CCMI data set is still used. Furthermore, also interactive aerosol activation into cloud droplets is calculated following the procedure described by . The activated aerosol number is subsequently used to determine prognostic cloud droplet number concentrations in the two moment cloud microphysics scheme of . Consequently, the whole large-scale cloud and condensation scheme of EMAC has been replaced by this alternative calculation. This has implications for cloud optical properties (first indirect aerosol effect), cloud lifetime (second indirect aerosol effect), rain and snow production (other indirect aerosol effects), and scavenging, and hence feedback on the chemical composition of the atmosphere. Overall, this simulation provides a more physical representation of the processes in the lower atmosphere at the expense of increased computational costs (see Table ).

Prescribed emissions

Anthropogenic emissions are incorporated as prescribed emission fluxes following the CCMI recommendations . Two data sets are considered: one data set is the MACCity

MACCity $=$ Monitoring atmospheric composition & climate (MACC)/City Zero (carbon) Energy (CityZEN)

emission inventory, which is applied for the RC1 simulations covering the period from 1950 to 2010. The second data set consists of a combination of ACCMIP and RCP 6.0 data . This data set is utilized for the long-term (hindcast and projection) simulations (RC2). In the RC1 simulations, emission data of the year 1960 have been repeatedly applied for the spin-up period (i.e. the simulated years 1950 to 1959); in RC2 the ACCMIP data have been applied from 1950 on.

Apart from the temporal coverage, the characteristic difference between both data sets is that MACCity considers a seasonal (monthly resolved) cycle, whereas ACCMIP and RCP 6.0 prescribe monthly values, which have been linearly interpolated from annual emission fluxes. Seasonality is only provided for biomass burning and ship emissions.

Emission data sets prepared for the simulations carried out here combine the broad range of sectors provided by the original underlying emission inventories into six categories, namely land, road, agricultural waste burning, shipping, aviation, and biomass burning. Moreover, the ground-based emission fluxes are distributed vertically to characteristic heights as described by . For MACCity the values of total NMVOCs (non-methane volatile organic compounds) for the anthropogenic sector were re-calculated from the corresponding species, since total NMVOCs were not provided by the original data set.

In addition to the anthropogenic emissions, some non-anthropogenic emissions such as NMHCs of biogenic origin, terrestrial dimethyl sulfide ($\mathrm{DMS}$), volcanic ${\mathrm{SO}}_{\mathrm{2}}$, ${\mathrm{NH}}_{\mathrm{3}}$, halocarbons, ${\mathrm{CH}}_{\mathrm{3}}\mathrm{I}$, and $\mathrm{OC}$ from secondary organic aerosol ($\mathrm{SOA}$, the latter only for the -aero- and -aecl-simulations) have been prescribed, mostly based on climatologies. Further details on the preprocessing steps for the emission data sets (e.g. speciation of total NMVOCs into individual species, description of the composition of the individual sectors, and the non-anthropogenic emissions) can be found in the document ESCiMo_emissions.pdf, which is part of the Supplement.

Annual total emissions of on-line calculated biogenic/soil $\mathrm{NO}$${}_x$ emissions in $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

In the RC1 simulations, prescribed emissions of the last available year (typically 2010 or 2011) have been used for later years as well, with one exception: in RC1SD-base-10a (see Table  and Sect. ) we used the RC2 set-up for the years 2012 and later.

On-line calculated emissions

Biogenic emissions

The emissions of $\mathrm{NO}$${}_x$ from soil and isoprene (${\mathrm{C}}_{\mathrm{5}}{\mathrm{H}}_{\mathrm{8}}$) from biogenic sources are calculated on-line using the submodel ONEMIS. Besides prescribed fields, like the distribution of cultivation and agriculture, the calculated $\mathrm{NO}$${}_x$ emissions depend mainly on the soil temperature and the soil wetness. The algorithm used is based on as described by . Estimates for the soil biogenic emissions of $\mathrm{NO}$ are highly uncertain. compared results of different models and satellite observation-based estimates. Most of the annual totals are in the range of 4 to 15 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$, not taking into account the uncertainties reported in the individual studies.

An overview of the annual totals for the different simulations is shown in Fig. . In all simulations the emission totals are between 5.4 and 6.3 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$ until the year 2010 and thus at the lower end of the estimated range from . While the simulations of the RC2 series are at the lower limit of this range, the simulations RC1SD-base-07 and RC1SD-base-08 (i.e. with specified dynamics including global mean temperature nudging) are at the top of this range, as they show a slightly higher soil temperature compared to the RC1SD-base-09/10. Until 2100 the RC2 simulations show an increase of the emissions up to about 6.8 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$. This trend is due to the increasing soil temperature, as the soil wetness does not show a trend and is essentially in dynamic equilibrium 2 years after initialization.

The simulated isoprene emissions depend on prescribed fields like the leaf area index and on-line calculated quantities from the base model (surface temperature and the net solar radiation). The algorithm is based on and implemented according to . For the biogenic isoprene emissions (see Fig. ) we simulate a range between 430 and 550 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$. Similar as for the biogenic $\mathrm{NO}$${}_x$ emissions the RC2 series is at the lower end and the RC1SD simulations at the upper end of this range. As the isoprene emissions are also dependent on the surface temperature, we see a similar increase of the emissions for the RC2 series up to about 650 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ until 2100. However, do not consider the ${\mathrm{CO}}_{\mathrm{2}}$ inhibition effect on isoprene emissions, implying that the future isoprene emissions (as projected in the RC2 simulations) may be too strong.

Annual total emissions of on-line calculated biogenic isoprene emissions in $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

estimated the annual total emissions of isoprene from biogenic origin to be 440 to 660 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ see alsoand references therein. All our simulations are within this range, however, these total emissions are further scaled with a factor of 0.6 to yield realistic mixing ratios of isoprene in the boundary layer see. A more detailed discussion on the isoprene scaling factor is provided by .

Lightning $\mathrm{NO}$${}_x$

The $\mathrm{NO}$${}_x$ emissions from lightning activity are calculated on-line using the MESSy submodel LNOX . Here, we apply the parameterization by , which links the flash frequency to the updraft velocity. The flash frequency obtained by this parameterization is scaled with 3.81459 for all simulations at L90MA and with 6.548 for all simulations at L47MA vertical resolution. As these scaling factors are applied identically for all simulation sets using the same vertical resolution, the total $\mathrm{NO}$${}_x$ emissions by lightning differ between the simulations.

Estimates for the annual total emissions of $\mathrm{NO}$${}_x$ from lightning are in the range of 2–8 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$ see. In Fig.  we see that most of the simulations are within this range. The results cluster around values between 4 and 5 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$. In the simulations with nudged global mean temperature (RC1SD-base-07/08), however, the total emissions are very low. The reason for this is a significantly reduced number of convective events due to a more stable temperature profile between the surface and the tropopause (see Sect. ). Despite the large difference of the annual total lightning $\mathrm{NO}$${}_x$ emissions between the RC1SD-base-10 and the RC1SD-base-07 simulations, the corresponding geographical distributions remain similar (see Figs. S25 and S26).

Annual total emissions of lightning $\mathrm{NO}$${}_x$ in $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

Ocean-to-atmosphere fluxes

Ocean-to-atmosphere fluxes of $\mathrm{DMS}$, ${\mathrm{C}}_{\mathrm{5}}{\mathrm{H}}_{\mathrm{8}}$, and methanol (${\mathrm{CH}}_{\mathrm{3}}\mathrm{OH}$) are calculated by the AIRSEA submodel . The ocean-to-atmosphere flux of a chemical species is calculated from its concentrations in the uppermost ocean layer and the lowermost atmosphere layer following the two-layer model by . We use the parameterization of for the water side exchange velocity in this study. The parameters needed for the simulated species in AIRSEA correspond to the values suggested by . Ocean salinity is taken from the monthly climatology of the World Ocean Atlas 2001 . Oceanic concentrations of $\mathrm{DMS}$ are a monthly climatology from . Isoprene concentrations in the ocean are calculated using the parameterization of , relating isoprene and chlorophyll concentrations, here with chlorophyll prescribed as a monthly climatology from the World Ocean Atlas 2001 . For methanol, the atmosphere-to-ocean flux is calculated assuming a constant undersaturation of surface water methanol with respect to lowermost atmospheric concentrations . The constant saturation coefficient for methanol in surface water is 0.94.

The annual global ocean-to-atmosphere fluxes of $\mathrm{DMS}$ are shown in Fig. ; total oceanic emissions range between 25.5 and 31.5 $\mathrm{Tg}\left(\mathrm{S}\right){\mathrm{a}}^{-\mathrm{1}}$ until 2010, and up to 32.5 $\mathrm{Tg}\left(\mathrm{S}\right){\mathrm{a}}^{-\mathrm{1}}$ in the projection simulations. These values are within the range reported in the literature of 15 to 54 $\mathrm{Tg}\left(\mathrm{S}\right){\mathrm{a}}^{-\mathrm{1}}$ and references therein. The trend in $\mathrm{DMS}$ emissions from the ocean follows the trend of SSTs (Fig. ). Therefore, compared to the -base- simulations, the increase in emissions is lower in the coupled ocean simulation RC2-oce-01, as the latter shows a lower increase of SSTs.

Simulated annual total emissions of $\mathrm{DMS}$ (in $\mathrm{Tg}\left(\mathrm{S}\right){\mathrm{a}}^{-\mathrm{1}}$) from the ocean.

[Figure omitted. See PDF]

Oceanic emissions of isoprene are considered low compared to biogenic emissions over land, mainly from tropical rainforests. Recent estimates for oceanic isoprene emissions range from 0.089 to 0.11 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ , much lower than previous estimates . Besides their low contribution to global totals, oceanic isoprene emissions are the main local isoprene source, in particular in the remote ocean marine boundary layer. Calculated oceanic isoprene emissions range from 0.063 to 0.074 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ (Fig. ), which is below the estimates found in literature, although still within the range of uncertainties. This underestimation is mostly due to the different chlorophyll distributions used to estimate the isoprene concentration in the surface water. The trend of isoprene emissions follows the trend of SSTs (Fig. ).

Simulated annual total emissions of isoprene (in $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$) from the ocean.

[Figure omitted. See PDF]

Simulated annual total flux of ${\mathrm{CH}}_{\mathrm{3}}\mathrm{OH}$ (in $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$) between atmosphere and ocean. The negative sign indicates a net uptake by the ocean.

[Figure omitted. See PDF]

For methanol, the net oceanic sink is uncertain, ranging from 0.1 to 21 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ , with a recent study by resulting in an net oceanic methanol sink of 6 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$. The simulated oceanic uptake of methanol ranges from 1.2 to 1.85 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ (Fig. ). These values are at the lower end of the range reported in earlier studies, which mostly assumed a fixed undersaturation of oceanic seawater of 90 $\mathrm{\%}$ , whereas we assumed 94 $\mathrm{\%}$. Furthermore, the net oceanic sink for methanol is still uncertain and only one publication reporting measurements of methanol concentrations in seawater exists .

Simulated annual total emissions of $\mathrm{Br}$ from sea salt (in $\mathrm{Tg}\left(\mathrm{Br}\right){\mathrm{a}}^{-\mathrm{1}}$).

[Figure omitted. See PDF]

Differences in the ocean-to-atmosphere fluxes between the free-running experiments and the experiments conducted with specified dynamics are mainly caused by differences in wind. The parameterization by relates the water side exchange velocity to the squared 10 m wind speed. Further differences occur because of different atmospheric composition and different atmospheric states between the experiments. The differences in the ocean-to-atmosphere methanol flux between the -aero-/-aecl- experiments and the -base- simulations are mainly caused by different atmospheric methanol concentrations, due to the applied constant undersaturation of surface water methanol.

We also simulate the release of $\mathrm{Br}$ from sea salt with submodel ONEMIS by scaling the mass flux of sea salt (accumulation and coarse mode) with the fraction of bromide in sea salt. Following , we assume that 50 $\mathrm{\%}$ of the bromide is released to the gas phase (leading to an additional factor of 0.5). The resulting time series of oceanic $\mathrm{Br}$ emissions is shown in Fig. . The absence of a trend in sea-salt emissions indicates that the 10 $\mathrm{m}$ wind speed over the ocean does not change significantly on global and annual average. This is consistent with the study by , revealing that changes in wind speeds from general circulation model (GCM) projections are small ($\le 15$ $\mathrm{\%}$).

Based on ocean organic carbon content (derived from SeaWiFS satellite data; see Fig. E57) an emission flux of oceanic particulate organic carbon (POC) associated with sea salt emissions in the accumulation mode is calculated in the -aero- and -aecl- simulations. This calculation (in submodel ONEMIS) follows the parameterization described by . The result is shown in Fig. . The oceanic organic carbon source is highly uncertain, with recent estimates of 8 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ and 75 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ . Our simulated annual total oceanic emissions of $\mathrm{POC}$ range between 17.9 and 19.5 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$, which is close to the first estimate of 14 $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$ for oceanic $\mathrm{POC}$ emissions by .

Simulated annual total emissions of $\mathrm{POC}$ (in $\mathrm{Tg}\left(\mathrm{C}\right){\mathrm{a}}^{-\mathrm{1}}$) from the ocean.

[Figure omitted. See PDF]

Dust emissions

Dust emissions are calculated on-line in the RC1-aero and RC1-aecl simulations. We use the dust emission scheme of see Figs. E58 and E59 for input parameters;, which was implemented in the submodel ONEMIS by . The wind stress threshold for dust emissions is corrected by a factor of 0.86 in accordance with . The on-line calculated total global dust emissions are in the range between 600 and 1400 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$ (Fig. ), of which 1.5 to 2 $\mathrm{\%}$ of the mass is emitted into the Aitken mode and the rest into the coarse mode. Whereas the dust emissions in the RC1-aero simulations are at the lower limit of the suggested range of 800 to 1700 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$ , dust emissions of the RC1-aecl simulations are within the suggested range.

As shown in Fig. , RC1-aero and RC1-aecl result in different dust emissions, which are sensible to wind speed, but also to surface dryness as a consequence of precipitation. The aerosol–cloud interactions modify the wind speed via boundary layer processes, which are induced by the differential heating caused by aerosol impacts on clouds. Additionally, the circulation is slightly altered, such that higher mean wind speed close to the surface is obtained. Additionally, precipitation (see Fig. ) is slightly different in the RC1-aecl simulation compared to the -base- case. For instance in Central Africa RC1-aecl is too wet compared to data from the Global Precipitation Climatology Project (GPCP), whereas the -base- case is underestimating precipitation slightly.

Annual total emissions of dust (in $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$) as calculated in the RC1-aero and RC1-aecl simulations.

[Figure omitted. See PDF]

The simulations with prescribed aerosol use the Tanre climatology (Sect. ), which explicitly accounts for mineral dust as one of the main components. Therefore, only the spectral climatological distribution of dust particles is used instead of emission fluxes.

Simulation with coupled ocean model

Simulation RC2-oce-01 with atmospheric chemistry and an interactively coupled ocean model covers the period 1950–2100. Based on the model set-up of the RC2-base-08 simulation, the submodel MPIOM was additionally switched on, together with the submodels A2O and HD. The dynamical coupling between ocean and atmosphere via the A2O submodel was computed every 2 h.

The simulation RC2-oce-01 with the coupled atmosphere–ocean EMAC set-up was performed at T42L47MA resolution for the atmosphere and at GR30L40 resolution for the ocean. This ocean model resolution corresponds to an average horizontal grid spacing of 3${}^{\circ }$ $\times$ 3${}^{\circ }$, with 40 unevenly spaced vertical levels. The rotated ocean model grid is shown in Fig. S6.

For the RC2-oce-01 simulation and the prior spin-up procedure (see below), we applied the parameter set as optimized by . In addition to the parameters described in Sect. , cloud optical property-related parameters have been modified (see also Table S2): the asymmetry factor of ice clouds $zasic=0.91$, the inhomogeneity factor of ice clouds $zinhomi=0.80$, and the inhomogeneity factor of liquid clouds $zinhoml=0.70$.

First, a sequence of two spin-up simulations has been conducted to provide internally consistent initial conditions of both, the ocean and the atmosphere component of the coupled system, representative for the year 1950. In these simulations only the dynamical components of EMAC-MPIOM were used, i.e. GHG mixing ratios (${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{CH}}_{\mathrm{4}}$, ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$, $\mathrm{CFCs}$) were prescribed according to the CCMI recommendation, and no interactive atmospheric chemistry was calculated. The first spin-up simulation (SP-oce-01) has been integrated over 300 $\mathrm{years}$ to reach a thermal and radiative equilibrated state between ocean and atmosphere for pre-industrial conditions with a reasonable radiation budget at the TOA. For simulation SP-oce-01, required to equilibrate the system for pre-industrial conditions, GHGs were prescribed representative for the year 1750, i.e. without any trend. The applied values are $0.28\times {10}^{-\mathrm{3}}$, $0.72\times {10}^{-\mathrm{6}}$, and $0.27\times {10}^{-\mathrm{6}}$ $\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}$ for ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{CH}}_{\mathrm{4}}$, and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$, respectively. The mixing ratios of the $\mathrm{CFCs}$ were set to zero. The remaining radiation imbalance in SP-oce-01 averaged over the last 30 $\mathrm{years}$ is $1.375$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ with a globally averaged surface temperature of 288.4 $\mathrm{K}$. This imbalance seems to be quite high; however, the corresponding optimized parameter set is the result of a multitude of sensitivity studies , yielding the optimum results in terms of climate and hydrological cycle for both, pre-industrial and industrial conditions. At the end of the spin-up period of SP-oce-01, no statistically significant trends (90 $\mathrm{\%}$ confidence level) of the surface temperature over the last 60 $\mathrm{years}$ were detectable.

From this equilibrated simulation SP-oce-01, a second simulation SP-oce-02 was integrated from the year 1750 to the year 1950 (i.e. 200 $\mathrm{years}$), with increasing (annually resolved prescribed) GHGs as suggested by CCMI, which pass over seamlessly into those used from 1950 onwards (see Figs. E8–E10). For both spin-up simulations, SP-oce-01 and SP-oce-02, ${\mathrm{O}}_{\mathrm{3}}$ has been prescribed from a monthly climatology, which has been derived from the years 1962–1972 of a previous test simulation. showed with a middle atmosphere chemistry climate model that “there are only minor changes in simulated stratospheric temperature and ozone prior to the year 1960”. Therefore, prescribing an ozone climatology for the simulated spin-up period might only introduce a small systematic error, as recently reported by .

Finally, simulation RC2-oce-01 was started from the end of SP-oce-02 (i.e. January 1950) with all additional submodels for atmospheric chemistry and diagnostics as in RC2-base-08.

Figure  shows the simulated global average SSTs and SICs of RC2-oce-01 compared to the corresponding time series prescribed in the other simulations. The results from RC2-oce-01 show an increase in temperature towards the end of the 21st century (see also Fig. S16) with respect to the 1980s, as is comparably projected by other numerical models Fig. AI.4, which show an ensemble median (average) increase of 1.7 (2.1) $\mathrm{K}$ for the same scenario (RCP6.0). Nevertheless, this increase (about 1.5 $\mathrm{K}$, Fig. , upper panel) is below what has been projected by the HadGEM-ES model (about 2.5 $\mathrm{K}$), which also implies a slower decrease of the global sea-ice coverage (Fig. , lower panel). Despite this smaller increase in temperature, the results are still in line with the CMIP5 simulations , where the projected increase is between 1.3 and 2.7 $\mathrm{K}$. Moreover, the HadGEM model is at the high end of the CMIP5 models in terms of climate sensitivity .

Global annual average (ice-free) SSTs (in $\mathrm{K}$, upper panel) and SIC (in ${10}^{12}{\mathrm{m}}^{\mathrm{2}}$, lower panel) of simulation RC2-oce-01 (1950–2099, black) compared to prescribed SSTs/SICs of simulations RC1-base (HadISST, blue), RC2-base (HadGEM2-ES, red) as described in Sect. , and of ERA-Interim (green) as applied in the nudged simulations RC1SD-base.

[Figure omitted. See PDF]

Glitches and unintended sensitivity studies

Due to the complexity of the model and the various model set-ups, gremlins have had a large zone of attack to creep in . Four issues have been detected during the course of the simulations or right after they had been finished. Given the large demand of resources (see Table ) and the advanced stage of the project, we were not able to repeat all simulations. Moreover, we had to ponder between achieving time series, which are as consistent as possible between the different simulations (despite their shortcomings), or to end up with “broken” time series. Except for the OC/BC emission issue detected in the RC1-aero and RC1-aecl simulations (Sect. ), we decided for consistent time series, but performed additional sensitivity simulations (indicated by suffix a; see Table ) to estimate the effects. These issues, described in more detail below, need to be taken into account in future analyses of the data.

Stratospheric aerosol optical properties

Due to a unit conversion error at data import, the extinction of stratospheric aerosols (Sect. ) was too low, by a factor of approximately 500. Since the contribution of stratospheric background aerosol to the radiative heating rate is minor, this is not a big issue. It has been tested by calculating the multi-annual monthly average radiative heating rates (1990–2010) from simulations RC1-base-08, RC1-base-08a, RC1-base-07, and RC1-base-07a. Above 100 $\mathrm{hPa}$ the range of differences between corresponding pairs is smaller than the inter-annual standard deviation. But very important: the dynamical effects of large volcanic eruptions (e.g. Mt. Pinatubo 1991; El Chichón 1982) are essentially not represented in the simulations, except for the contribution to the tropospheric temperature signal induced by the prescribed SSTs. The effect of stratospheric volcanic aerosol on infrared radiative heating is weak, as shown by mostly insignificant differences between RC-base-07a and RC-base-07, and RC-base-08a and RC-base-08, respectively (not shown).

The chemical effects (through heterogeneous chemistry, Sect. ), however, are included, since the prescribed aerosol surface areas were treated correctly. This gives rise to a very specific sensitivity study, which will be analysed elsewhere.

Road traffic emissions

Due to a wrong namelist entry, the timing of the road traffic emissions was unfortunately wrong: instead of updating the monthly input fields every month, they have been updated only every year, thus in 1950 emissions of January 1950 have been used, in 1951 the emissions of February 1950, etc. The resulting wrong emission time series are displayed in Figs. E19–E25 and E48–E54. All simulations are affected, except for RC1SD-base-10a, which constitutes a sensitivity simulation with both this and the stratospheric aerosol optical properties issue (Sect. ) corrected. The corrections cause an increase of the total ozone column of up to 2 $\mathrm{DU}$, mostly in the troposphere (see Fig. S17). It is expected that the effect on ozone is largest close to the road traffic emission regions. This has to be taken into account in future analyses.

Emissions of black carbon and organic carbon

Due to an incorrect unit conversion, the aerosol emissions for BC and OC in RC1-aero-06 and RC1-aecl-01 are substantially underestimated in terms of mass. However, particle numbers are correctly emitted in the corresponding number tracers. This results in a substantial underestimation of the aerosol radius, especially in the Aitken mode, and hence aerosol extinction, absorption, and activation are also too small in these simulations. As the feedbacks on chemically reactive compounds (both gas phase and secondary inorganic aerosol species) are only via the diffusion limitation, the consequences for other chemical species are minor. The fine mode aerosol distributions are, however, quite substantially impacted by the errors in OC/BC. Since the total budgets of many compounds are dominated by the larger size categories, they, except for OC/BC, are hardly affected. Furthermore, the impact on the aerosol optical properties of the small particles is also lower than for larger particles, such that the impact on radiation is also minor. For aerosol–cloud interactions (-aecl-) the error is only in the very first phase of the simulation leading to an underestimation of cloud droplets. As the problem has been fixed before the dominant change in especially organic aerosol emissions, the effect of increased cloud droplets from the year 1970 onwards is included in the resulting time series of RC1-aecl-02. As a consequence, detailed analyses of OC/BC can safely be based on results from RC1-aero-07 (from 1991 onwards) and RC1-aecl-02 (from 1966 onwards), respectively.

Heterogeneous ice nucleation

RC1-aecl-01 and RC1-aecl-02 overestimate the heterogeneous ice nucleation, due to an error in the calculation of ice nucleus numbers. Even though this has implications for the radiation budget influenced by cirrus clouds, the effect is minor (as shown by sensitivity tests after correcting the error).

Selected results

This section presents some first analyses, intended only to serve as a first overall evaluation and consistency check of the simulation results. Note that we apply in our analyses observational data without considering uncertainties in measurements, retrievals or uncertainties arising from representativity e.g. and without a detailed analysis of statistical implications on the evaluation of the simulation data e.g.. As often as possible, we added information on inter-annual variability or sampled data exactly as it is done for the observational data in order to reduce some of the uncertainties in the evaluation. More in-depth analyses of specific topics will follow and will be published elsewhere.

Upper left panel: climatology of annual average total dry air temperature of ERA-Interim in $\mathrm{K}$. The data were monthly and zonally averaged for the period 2000–2010. Other panels: dry air temperature differences (in $\mathrm{K}$) of the simulations compared to ERA-Interim data. The differences, unless grey shaded, are significant on a 95 $\mathrm{\%}$ confidence level according to a two-sided Welch's test.

[Figure omitted. See PDF]

Temperature profiles

For the assessment of the simulated air temperatures, we compare the simulation results with ERA-Interim data . We choose the period 2000–2010 for the evaluation, as it is covered by most of the simulations. Multi-annual climatologies of zonally averaged temperature profiles for the period 2000–2010 are shown in Fig. . Figures  and S18 were produced with the ESMValTool . The data were first monthly then annually averaged and the values of ERA-Interim were subtracted from the simulation results. For the discussion, the simulations with basically identical characteristics concerning the temperature, can be put into three categories: First, the nudged simulations including nudged global mean temperature: RC1SD-base-07 and RC1SD-base-08. Second, the nudged simulations (without global mean temperature nudging): RC1SD-base-09 and RC1SD-base-10. Third, the free-running simulations: RC1-base-07, RC1-aecl-02, RC2-oce-01, and RC2-base-04. Here, we only discuss the temperature biases, the effect on stratospheric water vapour is investigated elsewhere .

The simulations in the first category perform best, as expected, compared to ERA-Interim with minor differences of mostly less than $\pm$1 $\mathrm{K}$. Furthermore, there are only minor differences visible comparing the two different vertical resolutions of the simulations. For this category we only found a small seasonal variation (Fig. S18). Solely, in the Northern Hemisphere (NH) winter the simulations show a particular warm bias in the north polar region.

The simulations in the second category show a cold bias, which moreover has an obvious vertical structure. However, the vertical patterns of the simulations in this category differ significantly. The maximum of the cold bias in the RC1SD-base-09 simulation is at around 70 $\mathrm{hPa}$ with more than $\pm$4 $\mathrm{K}$. Between 100 and 150 $\mathrm{hPa}$ this simulation shows a small bias of $\pm$1 to $\pm$1.5 $\mathrm{K}$. In the RC1SD-base-10 simulation, with a higher vertical resolution, the maximum of the cold bias is located lower, at about 200 $\mathrm{hPa}$. The differences in the patterns are an effect of the vertical resolution, but – although significant – the nudging of the model has a much larger impact on the temperature distribution (cf. RC1SD-base-10 with RC1-base-07 in Fig. ). Concerning the seasonal variations, all simulations in this category show a similar behaviour, as the temperature bias in December, January, and February is the smallest in the Arctic, while it is overall mostly negative.

The simulations in the third category have a cold bias (of up to 4.5 $\mathrm{K}$) around the tropopause and show a warm bias (of up to 4 $\mathrm{K}$) in the Southern Hemisphere (SH) above 100 $\mathrm{hPa}$, as also analysed from previous EMAC simulations by .

In this region the bias is smallest in the RC2-oce-01 simulation, but the data show a larger cold bias below the tropical (30${}^{\circ }$ S–30${}^{\circ }$ N) tropopause at around 250 $\mathrm{hPa}$. The free-running simulations show a more pronounced seasonal cycle of the bias than the nudged simulations. The warm bias above the tropopause around 60–30${}^{\circ }$ S, which is also seen in the annual climatology, is strongest in SH winter (more than 4 $\mathrm{K}$). Around this time there is also a strong cold bias (of up to 4 $\mathrm{K}$) in the tropics above the tropopause.

(a) Zonal distribution of mean precipitation ($\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) for the time period 1990–2009 for all simulations covering this period plus GPCP and ERA-Interim. (b) Deviation of the simulated zonal distribution of precipitation ($\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) shown as simulation results minus GPCP and ERA-Interim minus GPCP. Colour-coding as in panel (a). Shading in both panels (a) and (b) is $\pm$ 1 $\mathit{\sigma }$ of the annual variation in GPCP data (1990–2009). The numbers in panel (b) list the corresponding global and tropical (30${}^{\circ }$ S–30${}^{\circ }$ N) average precipitation (in $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$), respectively.

[Figure omitted. See PDF]

Hydrological cycle

To evaluate the results of the simulations with respect to the representation of the hydrological cycle, we compare them with global precipitation data from the GPCP. Additionally, we include total precipitation from ERA-Interim reanalysis http://apps.ecmwf.int/datasets/data/interim-full-mnth/; to bridge the gap between models and observations and to have a reference for the specified dynamics simulations. GPCP version 2.2 http://www.esrl.noaa.gov/psd/data/gridded/data.gpcp.html; data are available as monthly means for each month since January 1979 at 2.5${}^{\circ }$ $\times$ 2.5${}^{\circ }$ horizontal resolution. This precipitation data set combines satellite data with rain gauge measurements, resulting in a data set with global coverage. As , we have used cumulated precipitation at step 0 and step 12 of the reanalysis' total precipitation. For the simulations we used the “snapshots” of the precipitation fluxes (large-scale and convective rain and snow fall) to calculate the flux of total precipitation. The period 1990–2009 was chosen in order to be able to analyse as many of the simulations as possible. First, we look at zonally averaged precipitation in the years 1990–2009 as shown in Fig. .

(a–f) Mean precipitation differences ($\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) for the 20-year period 1990–2009. The differences show the simulation results minus GPCP data.

[Figure omitted. See PDF]

All simulations, as well as ERA-Interim, overestimate precipitation in the tropics by 0.14 up to 0.84 $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$, depending on the simulation (see Fig. ). However, simulations and reanalysis show a double-peak structure (with varying intensity) near the Equator with a stronger peak (6 to 7.5 compared to 4 to 5 $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) in the NH and minima (about 2 $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) of precipitation at approximately 25${}^{\circ }$ N and 25${}^{\circ }$ S as is also apparent in the GPCP data. The maxima in the SH and NH storm tracks as shown by GPCP (about 3 and 2.75 $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) are reproduced by ERA-Interim and all EMAC simulations (about 3 to 3.75 and 2.25 to 2.75 $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$). Apart from the tropics, the largest deviations can be found in the SH subtropics. Horizontal distributions of average precipitation deviations (in $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) throughout the years 1990–2009 are shown in Fig. .

We show deviations as model results minus GPCP exemplarily for six simulations (SD with global mean temperature nudging, SD without mean temperature nudging, free running, free running with coupled aerosol, transient, and transient with coupled ocean). Simulation pairs which differ only in vertical resolution show very similar precipitation patterns. Moreover, precipitation patterns and zonal distribution of precipitation from the simulations with corrected stratospheric aerosol (RC1-base-07a/08a; brown lines in Fig. ) and the simulation with tropospheric aerosol (RC-aero-07, orange line) are very similar to the respective base cases (RC1-base-07/08 and RC1-base-07, respectively; red lines).

Corresponding results of all other simulations and the GPCP data are shown in Figs. S19 and S20. Nudged simulations (RC1SD-base; see also cadet blue and blue lines in Fig. ) show slightly lower deviations from observations than free-running simulations. In general, the free running and the nudged simulations show the same large-scale deviation patterns from the observations. However, they differ in strength and also regional differences can be found. The largest absolute deviations (of more than $\pm$3 $\mathrm{mm}{\mathrm{day}}^{-\mathrm{1}}$) are found in the tropics, where rainfall is the strongest.

Globally integrated, annual dry deposition flux of ozone in $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

Overall the magnitude of deviations is comparable to uncertainties that arise when changing convection parameterizations within ECHAM5 and have also been found in the analysis of the ECHAM5 model . The simulation which differs the most is RC2-oce-01 (dashed purple line), i.e. the simulation with coupled ocean. In particular in the tropics the simulation produces a double Inter Tropical Convergence Zone (ITCZ), a typical feature of coupled models .

Deposition fluxes

In EMAC three processes establish the sink of trace gases and aerosols to the surface: dry and wet deposition and (for aerosols only) sedimentation. As an example for gas-phase species, Fig.  displays the dry deposition flux of ozone. Generally, all simulations show the same temporal evolution. Starting in 1950 with deposition fluxes between approximately 700 to 730 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$, the deposition fluxes increase up to a maximum of 850 to 880 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$ around the year 2000. This is well within the range of other chemistry climate model results. For instance, reported ozone deposition fluxes in the range between 687 and 1360 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$ for six models taking part in the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP). The dry deposition fluxes of the RC1-base and RC1SD-base simulations stay in this range until the end of the simulations. In the projection simulations (RC2) the ozone dry deposition fluxes decrease from the year 2060 onward. The temporal evolution of the deposition fluxes mirrors the temporal evolution of the tropospheric ozone mixing ratio (not shown). The ozone dry deposition fluxes are higher in the simulations with lower vertical resolution (L47MA) compared to L90MA. Striking are the two RC1SD simulations with “wave zero” temperature nudging. They produce much lower ozone dry deposition fluxes, which is a direct effect of the largely reduced ozone mixing ratio (about 28 to 32 $\mathrm{nmol}{\mathrm{mol}}^{-\mathrm{1}}$ on average) in the RC1SD-base-07/08 simulations compared to other simulations. The lower ozone mixing ratio, in turn, is caused mainly by the reduced lightning $\mathrm{NO}$${}_x$ and corresponding ozone production (see Fig.  and Sect. ). The wet deposition fluxes of ozone exhibit the same features, but are 6 to 7 orders of magnitude smaller (see Fig. S21).

For the $\mathrm{NO}$${}_x$ wet deposition fluxes the difference between L47MA and L90MA amounts to 20 $\mathrm{\%}$ (see Fig. S22). In contrast to ozone, the $\mathrm{NO}$${}_x$ dry deposition fluxes peak earlier (around 1990) and decrease slowly afterwards (see Fig. S22).

The deposition fluxes for nitrate (${\mathrm{HNO}}_{\mathrm{3}}+{\mathrm{N}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{5}}+{\mathrm{NO}}_{\mathrm{3}}^{-}$) and sulfate (${\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}+{\mathrm{HSO}}_{\mathrm{4}}^{-}$ $+$ ${\mathrm{SO}}_{\mathrm{4}}^{\mathrm{2}-}$) are shown in Figs.  and , respectively. Fluxes of ammonium and sulfite (species with sulfur in oxidation state $+$4, $\mathrm{S}\left(\mathrm{IV}\right)$) are shown in Figs. S23 and S24, respectively. All RC1-base and RC1SD-base simulations behave similarly. For nitrate all deposition fluxes increase between 1950 and 2010, while the sulfate deposition fluxes peak somewhere between 1970 and 2000. The deposition fluxes of the projection simulations are lower in the 50s. For the dry deposition and the sedimentation fluxes the vertical resolution seems to make a noticeable difference, while for the wet deposition the presence of an interactive ocean has a larger impact on this sink. However, the dominant sink for nitrate and non-sea salt (nss) sulfate is wet deposition.

Globally integrated annual nitrate dry deposition (upper), wet deposition (middle), and sedimentation (lower) fluxes in $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

Globally integrated annual non-sea salt sulfate dry deposition (upper), wet deposition (middle), and sedimentation (lower) fluxes in $\mathrm{Tg}\left(\mathrm{S}\right){\mathrm{a}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

The largest differences can be found for the simulations including aerosol chemistry. The nitrate dry deposition fluxes and the sedimentation fluxes are larger, if aerosol chemistry is included. This results from the size of the aerosol particles. While the radii are prognostically calculated in the simulations with aerosol chemistry, they are simply prescribed in the other simulations. As the simulated coarse-mode radii are greater than the prescribed ones (globally averaged coarse-mode aerosol radii: 0.5 $\mathrm{\mu }$m (prescribed); simulated: 0.759 $\mathrm{\mu }$m (RC1-aero-06), 0.751 $\mathrm{\mu }$m (RC1-aero-07), 0.793 $\mathrm{\mu }$m (RC1-aecl-01), 0.782 $\mathrm{\mu }$m (RC1-aecl-02)), the dry deposition and sedimentation velocities are larger. Nevertheless, the wet deposition fluxes are smaller if aerosol chemistry is taken into account, yielding – summed over all deposition processes – the same amount of deposited nitrate.

Although the radii in the simulations with additional consideration of aerosol cloud interactions (-aecl-) are the largest, the deposition fluxes for all three deposition processes are slightly smaller. The largest effect is found on aerosol sedimentation as the sedimentation velocity depends the strongest on the particle radius .

For the sulfate deposition fluxes, the differences between the -base- and the -aero-/-aecl- simulations are much larger than for nitrate. This originates in a different treatment of sulfate emissions in the simulations with aerosol chemistry. In these simulations sea salt sulfate is taken into account as additional sulfate source. To keep the sulfate depositions comparable, only nss-sulfate is depicted in Fig.  and discussed below. The additional emissions, of course, influence the properties of the aerosol particles and thus differences had to be expected. Again, the nss-sulfate sedimentation fluxes (Fig. , lower panel) noticeably depend on the vertical resolution. The dry deposition of all -base- simulations are very similar except for the first 10 years, when the projection simulations are distinctively lower than the other simulations. The largest differences occur for the simulations regarding aerosol chemistry. The dry deposition fluxes are up to 50 $\mathrm{\%}$ smaller, the wet deposition fluxes up to 15 $\mathrm{\%}$ smaller in -aero- and -aecl- simulations, while the sedimentation flux is up to 2 times larger. The overall sink of nss-sulfate due to the three deposition processes is approximately 10 $\mathrm{\%}$ smaller in the simulations with aerosol chemistry.

The deposition fluxes for sulfur and oxidized and reduced nitrogen approximately balance the emissions for present-day values. Compared to the ACCMIP evaluation for the year 2000 sulfur emissions (especially $\mathrm{DMS}$) are slightly enhanced, whereas $\mathrm{NO}$${}_x$ and ${\mathrm{NH}}_{\mathrm{3}}$ emissions are a bit lower. The ratio of dry deposition to wet deposition is slightly lower in the present simulations compared to the model mean determined by , with, e.g., 0.64 for oxidized sulfur compared to 0.72. Similar differences are found for $\mathrm{NH}$${}_x$ and $\mathrm{NO}$${}_{x,y}$, with ratios of 0.5 compared to 0.66 and 0.63 to 0.72, respectively. This is most likely due to the more explicit treatment of wet deposition and cloud and precipitation chemistry in the EMAC model. Previous model versions have shown a reasonable agreement with surface deposition in North America, Europe, East Asia, and Africa .

Tropospheric oxidation capacity

The hydroxyl radical ($\mathrm{OH}$) is the most important cleansing oxidant in the atmosphere. However, variations in $\mathrm{OH}$ concentration are still difficult to quantify . Among others, it reacts with ${\mathrm{CH}}_{\mathrm{4}}$ largely determining its atmospheric lifetime. Here, we present the simulated $\mathrm{OH}$ lifetime of atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ at time $t$ as a measure for the oxidising power of the atmosphere, calculated according to as $$\begin{array}{cc}& {\mathit{\tau }}_{{\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}}\left(t\right)=\\ & & {\displaystyle \frac{\sum _{b\in B}{M}_{{\mathrm{CH}}_{\mathrm{4}}}^b\left(t\right)}{\sum _{b\in B}{k}_{{\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}}^b\left(t\right)\cdot c_{\text{air}}^b\left(t\right)\cdot {\mathrm{OH}}^b\left(t\right)\cdot M_{{\mathrm{CH}}_{\mathrm{4}}}^b\left(t\right)}},\end{array}$$ where $M_{{\mathrm{CH}}_{\mathrm{4}}}^b$ being the mass of methane, $k_{{\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}}^b$ the reaction rate of the reaction ${\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}$, $c_{\text{air}}^b\left(t\right)$ the concentration of air, and ${\mathrm{OH}}^b\left(t\right)$ the mole fraction of $\mathrm{OH}$ in the grid box $b\in B$, with $B$ being the set of all grid boxes, which are below the on-line calculated tropopause. The lifetime of atmospheric methane is determined on the one hand by the abundance of $\mathrm{OH}$ and on the other hand by temperature, since the reaction rate of ${\mathrm{CH}}_{\mathrm{4}}+\mathrm{OH}$ depends on temperature. Figure  shows the simulated tropospheric methane lifetimes (towards $\mathrm{OH}$). It is first calculated for every output time step (i.e. 10-hourly), then averaged monthly, and then annually. Overall the simulated lifetime increases between 1950 and 1975 and decreases until 2013. The RC2-base-04 and RC2-base-05 simulations (orchid, solid, and dashed lines, respectively) predict a similar decreasing behaviour for the future. The simulation with the interactively coupled ocean model RC2-oce-01 (purple) behaves similarly, but predicts at first a small increase until 2070 before it decreases. The tropospheric methane lifetime of all simulations varies around an average value of $8.0\pm 0.6$ $\mathrm{years}$ in the period of 2000–2004 (mean and standard deviation for all simulations covering this period). In comparison to other methane lifetime (versus $\mathrm{OH}$) estimates, as by , , , , and ($9.8\pm 1.6$, $8.02$, $8.7$, $8.3$, $7.9$–$9.1$ years, respectively), the results of the current study tend to be lower, but mostly stay within the uncertainty range. Observation-based estimates derived from methyl chloroform abundance are generally longer than the model-based estimates e.g.with $11.2\pm 1.3$ and $10.2$ ($+$0.9/$-$0.7) years, respectively. The wide range of lifetime estimates is mainly caused by different methods of calculation and applied weighting , whereas varying included vertical layers due to different tropopause heights have a minor impact see also. The global mean temperature nudged RC1SD-base-08 simulation (royalblue dashed) predicts the shortest lifetime. The nudging leads to a higher temperature in basically all vertical layers compared to the simulations without global mean temperature nudging (see Fig. S18), which accelerates the oxidation of ${\mathrm{CH}}_{\mathrm{4}}$ by $\mathrm{OH}$ and also has a positive impact on the $\mathrm{OH}$ production by tropospheric ozone. This, however, does not hold for RC1SD-base-07 (royalblue solid) compared to RC1SD-base-09 (cadet blue dashed). The lifetime in RC1SD-base-07 is longer than in RC1SD-base-09, although average temperatures are higher. This results from a lower $\mathrm{OH}$ mixing ratio in RC1SD-base-07, mainly in the tropical upper tropopause, caused by a likewise smaller $\mathrm{NO}$${}_x$ production from lightning activity (see Fig. ). The longest lifetime is simulated in RC1-aecl-02 (burlywood), which shows the largest cold biases of all simulations.

Annually averaged methane lifetime with respect to $\mathrm{OH}$ of the different simulations calculated with Eq. ().

[Figure omitted. See PDF]

Moreover, the variations in production of $\mathrm{OH}$ by tropospheric ozone depend on various factors . However, a detailed discussion of these variations and the impact on the methane lifetime is beyond the scope of this study.

Chemistry in the upper troposphere and lower stratosphere (UTLS): comparison with CARIBIC data

In this section, we present results comparing model output of the five nudged simulations (RC1SD-base) with measurements of the project IAGOS-CARIBIC (In-service Aircraft for a Global Observing System – Civil Aircraft for Regular Investigation of the Atmosphere Based on an Instrument Container). As described by , an air freight container equipped with a number of instruments is deployed on a civil aircraft on three to four intercontinental flights per month. The project has been ongoing in its second phase since the year 2005 and the nudged simulations cover the period until the end of 2013, so measurement data of 348 flights are available for comparison with the model results. The 10-hourly model output was interpolated linearly in latitude, longitude, logarithm of pressure and time to the location of the aircraft by post-processing in order to work with comparable data sets of measurement and model data.

We have chosen to compare the model results and measurements in the form of seasonal climatologies, using data of ozone, methane, carbon monoxide, and acetone ($({\mathrm{CH}}_{\mathrm{3}}{)}_{\mathrm{2}}\mathrm{CO}$). CARIBIC flies to destinations worldwide; thus in order to obtain meaningful climatologies with good data coverage, latitude was limited to 35${}^{\circ }$ N $<$ latitude $<$ 60${}^{\circ }$ N. Data were also limited by including only values where pressure $p<280$ $\mathrm{hPa}$ to exclude ascents and descents of the aircraft.

Climatologies are presented with a vertical scale relative to the tropopause in kilometres. This information is provided with the measurement data, and derived by linear interpolation to the aircraft position from ECMWF operational analysis data (with a 1${}^{\circ }$ grid resolution) as distance to the 3.5 $\mathrm{PVU}$ potential vorticity iso-surface (description available under http://www.knmi.nl/samenw/campaign_support/CARIBIC/). For comparability, i.e. to have the same tropopause definition, a similar calculation was applied to the model data. It is essential to reference data to this height relative to the tropopause (HrelTP), as stratospheric and tropospheric values then become separable when using the data from an aircraft that flies at constant pressure.

Comparison of ${\mathrm{O}}_{\mathrm{3}}$ climatologies (35–60${}^{\circ }$ N) based on data from the years 2005–2013. The first row shows climatologies of CARIBIC measurements (left) and RC1SD-base-10a model data (right). The second row shows the relative differences of model data and measurements (left), and the absolute values of model minus measurements over the sum of the standard deviation of measurements and model (right). The bottom row shows the relative differences of RC1SD-base-07 and RC1SD-base-10, differing in nudging of global mean temperatures (left), and the relative differences of RC1SD-base-07 and RC1SD-base-08, with different vertical resolutions (right). The vertical axis shows the distance (in $\mathrm{km}$) relative to the tropopause.

[Figure omitted. See PDF]

As Fig. , but for $\mathrm{CO}$.

[Figure omitted. See PDF]

As Fig. , but for ${\mathrm{CH}}_{\mathrm{4}}$.

[Figure omitted. See PDF]

As Fig. , but for acetone. The maximum of measured acetone values in summer is 1.5 $\mathrm{nmol}{\mathrm{mol}}^{-\mathrm{1}}$.

[Figure omitted. See PDF]

Results of the data prepared in this way are presented in Fig.  for ${\mathrm{O}}_{\mathrm{3}}$, Fig.  for $\mathrm{CO}$, Fig.  for ${\mathrm{CH}}_{\mathrm{4}}$, and Fig.  for acetone. The figures show the climatologies of measurements of RC1SD-base-10a and their relative differences. In order to judge upon the magnitude of the relative differences, a plot of the absolute values of model results minus measurements over the sum of the standard deviation of measurements and model results is also displayed. To investigate the different model simulations, two more climatologies are presented for each species: RC1SD-base-07 minus RC1SD-base-08, differing in vertical resolution, and RC1SD-base-07 minus RC1SD-base-10, differing in nudging of the global mean temperature. The results of RC1SD-base-10 minus RC1SD-base-10a, using different road emissions, are discussed, but shown only in the Supplement (Figs. S27–S34). In addition, the Supplement includes the climatologies and relevant relative differences of all nudged simulations (RC1SD-base). The subsequent paragraphs present the results separately for each species.

In the case of ${\mathrm{O}}_{\mathrm{3}}$ (Fig. ), the model results compare well with measurements (see , for instrument specifics). The seasonal cycle is reproduced, but the lower values in the troposphere are generally overestimated by up to 40 $\mathrm{\%}$ by the model. In the stratosphere, differences are smaller, as the model underestimates measurements by 5 $\mathrm{\%}$, reaching 30 $\mathrm{\%}$ only in summer (June through September), where the elevated ozone levels drop faster in the model than indicated by measurements. This earlier drop is especially noteworthy since it is larger in magnitude than the sum of the standard deviation of model and measurements. Correcting the road traffic emissions has only a minor influence on the distribution of ${\mathrm{O}}_{\mathrm{3}}$ (shown in Supplement). When nudging global mean temperature, upper tropospheric ${\mathrm{O}}_{\mathrm{3}}$ levels increase by about 5 $\mathrm{\%}$, while stratospheric levels stay constant. Increasing the vertical resolution on the other hand seems to decrease the differences between the seasons, reducing ${\mathrm{O}}_{\mathrm{3}}$ levels in spring and summer by about 10 $\mathrm{\%}$ in the simulations with 90 levels and increasing them by up to 5 $\mathrm{\%}$ in the second half of the year. When comparing results of this section with those of Sect. , it is important to keep in mind that the vertical scale is different and the selection of years is also not the same.

For $\mathrm{CO}$ (Fig. ; see , for instrument specifics), results mostly mirror those of ${\mathrm{O}}_{\mathrm{3}}$. The seasonal cycle is again reproduced by the model, showing a faster drop in the troposphere in spring, which leads to an underestimation of about 35 $\mathrm{\%}$, i.e. the model values deviate more than 1 standard deviation from the measurements. Stratospheric values are overestimated by about 10 $\mathrm{\%}$ throughout the year. Correcting the road traffic emissions has little influence on the results, values increase only by about 5 $\mathrm{\%}$ in the winter stratosphere (shown in Supplement). Nudging global mean temperature has a small effect on $\mathrm{CO}$: for the simulations with 90 levels (RC1SD-base-07 and 10), stratospheric springtime values (March to June) increase by about 8 $\mathrm{\%}$. Contrary to ${\mathrm{O}}_{\mathrm{3}}$, it is the increasing vertical resolution that produces differences between stratosphere and troposphere for $\mathrm{CO}$, tropospheric values increasing by 5 $\mathrm{\%}$ for the simulations with 90 levels, and being reduced by up to 10 $\mathrm{\%}$ in the stratosphere.

Results for ${\mathrm{CH}}_{\mathrm{4}}$ (Fig. ; see , for instrument specifics) are similar. Relative differences in general are very small for ${\mathrm{CH}}_{\mathrm{4}}$, reaching only 5 $\mathrm{\%}$, if model results are compared to measurements, and less than 1 $\mathrm{\%}$ among the different model simulations. The model results are in general too low compared to measurements, with a relative difference of $-$0.5 $\mathrm{\%}$ in the troposphere. The seasonal cycle is correctly reproduced, with slightly overestimated early summer stratospheric values (May to July) by about 2 $\mathrm{\%}$ and an underestimation of about 4 $\mathrm{\%}$ during the rest of the year; i.e. these differences are almost all larger than the sum of model and measurement standard deviation. As ${\mathrm{CH}}_{\mathrm{4}}$ is prescribed by Newtonian relaxation at the surface (see Sect. ), the values at the surface always correspond to measurements. This explains why the differences to measurements, even in the UTLS, are so small. Differences between the model simulations indicate how they differ in vertical transport or in methane lifetime (see Sect. ). ${\mathrm{CH}}_{\mathrm{4}}$ is noticeably affected by the road traffic emissions, decreasing the mean values by about 1 $\mathrm{\%}$ randomly, not showing any seasonality or stratification (shown in Supplement). This is consistent with a decreased ${\mathrm{CH}}_{\mathrm{4}}$ $\mathrm{OH}$ lifetime (Fig. ). In contrast, a vertical dependence becomes visible when taking the relative difference of simulations with different vertical resolution, where values are up to 0.2 $\mathrm{\%}$ in the troposphere and down to $-$0.5 $\mathrm{\%}$ in the stratosphere for model simulations with an increased vertical resolution (see Sect. ). Nudging of global mean temperature has an interesting effect on the results. Values are practically identical for all seasons, except for late winter (January to March) tropospheric and summer (March to October) stratospheric values, where nudging global mean temperature leads to a relative decrease of up to 0.8 $\mathrm{\%}$.

Acetone (Fig. ) shows the largest deviations, when comparing model results with measurements for instrument specifics, see. Tropospheric values are underestimated by about 60 to 100 $\mathrm{\%}$ and the seasonal cycle of the measurement data is not visible. The seasonal cycle is, however, reproduced when taking more model data from the UTLS into account, including data from longitudes different to those where CARIBIC flies (not shown). Again, the difference to measurements in the troposphere is larger than 1 standard deviation of model and measurements. The effect of different model simulations is small compared to this difference. Correcting the road traffic emissions has no influence on the results. Nudging global mean temperature increases mean values by about 15 $\mathrm{\%}$ in the stratosphere during late winter and spring (January to May). Increasing the vertical resolution also increases mean acetone values by 5 $\mathrm{\%}$ in all seasons and at all heights, except for summer and autumn (May to November) stratospheric values, which are reduced by up to 10 $\mathrm{\%}$. These results are in line with the previous evaluation of Table 8, who observed an underestimation of $\sim 50$–35 $\mathrm{\%}$ for the same biogenic emissions as used here (55 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$), which account for $\sim 85$ $\mathrm{\%}$ of the total emissions. These results indicate the need for considering oceanic emissions and photochemical production of ${\mathrm{CH}}_{\mathrm{3}}{\mathrm{COCH}}_{\mathrm{3}}$ from monoterpenes, methylbutenol, and higher iso-alkanes , which were not included. For example showed that $i$-butane and $i$-pentane are responsible for the photochemical production of 4.3 and 5.8 $\mathrm{Tg}{\mathrm{a}}^{-\mathrm{1}}$ of ${\mathrm{CH}}_{\mathrm{3}}{\mathrm{COCH}}_{\mathrm{3}}$, respectively.

By looking at the climatologies and their respective relative differences, it has become clear that there exist small but systematic differences between the results of the different model set-ups. The corrected road emissions have only a minor influence. Increasing the vertical resolution has effects that change mean values by up to 10 $\mathrm{\%}$. Different strengths of this influence can be noted for most species between the stratosphere and the troposphere, for ${\mathrm{O}}_{\mathrm{3}}$ also between different seasons. The effect of nudging global mean temperature is weaker, but it also shifts the seasonal cycle for all species except ${\mathrm{O}}_{\mathrm{3}}$, leading to relative differences of up to 10 $\mathrm{\%}$ for selected months. Overall, the model produces realistic values for ${\mathrm{O}}_{\mathrm{3}}$, but underestimates ${\mathrm{CH}}_{\mathrm{4}}$, $\mathrm{CO}$, and acetone, especially in the troposphere.

Stratospheric dynamics

This section deals with the stratospheric dynamics, i.e. how well the Brewer–Dobson circulation (BDC) is represented in the simulations. The BDC is the large-scale stratospheric and mesospheric transport circulation and determines amongst others the distribution of chemical trace gases. It refers to the large-scale residual circulation (RC) with upwelling in the tropics and sinking motion at mid- and high latitudes in the stratosphere and winter mesosphere. Moreover, quasi-horizontal mixing processes contribute to the stratospheric transport . The RC is often expressed in terms of the tropical upward mass flux ($F_{\text{trop}}$), which corresponds to the upward-directed mass transport in the tropics and is balanced by the downward mass fluxes in the extratropics of both hemispheres . The travel time of an air parcel from the entry region at the tropical tropopause to any region in the stratosphere is referred to as age of stratospheric air AoA;. It is obtained from linearly increasing surface emissions of an inert tracer, e.g. sulfur hexafluoride (${\mathrm{SF}}_{\mathrm{6}}$), which are mostly used from observations. The AoA contains both components of the BDC, RC and mixing processes.

Results on $F_{\text{trop}}$ and mean AoA are shown in the Sect.  and . For evaluation of the model results on $F_{\text{trop}}$, the RC is calculated from ERA-Interim reanalysis. Moreover, data of atmospheric ${\mathrm{SF}}_{\mathrm{6}}$ concentrations from July 2002 to April 2012 are available from the Michelson Interferometer for Passive Atmospheric Sounding MIPAS; instrument for the calculation of mean AoA. Here, we use mean AoA from the MIPAS-ENVISAT (Environmental Satellite) level 1b spectra version 5 , which is an update of .

The tropical upward mass flux

Figure  shows the time evolution of the annual mean tropical upward mass flux from 1960 to 2100 from the simulations in the lower (around 70 $\mathrm{hPa}$, top), middle (around 10 $\mathrm{hPa}$, middle), and upper (around 1 $\mathrm{hPa}$, bottom) stratosphere. ERA-Interim data are included from 1979 to 2013. All simulations show an increase in $F_{\text{trop}}$ with rising atmospheric mixing ratios of GHGs, which is well known from earlier CCM simulations e.g.. In the lower and middle stratosphere, the model simulations are grouped with respect to their underlying SSTs, whereas the vertical resolution has a minor influence on the performance of the RC. The RC2-simulations are in better agreement with ERA-Interim data than the simulations RC1-base-07/08 for the 1980s and 1990s. After 1995 the ERA-Interim data show a negative trend in the tropical upward mass flux (at 70/74 $\mathrm{hPa}$), which is not captured by the model simulations. As shown by , this negative trend contains large uncertainties, and does not occur in other reanalysis systems or when using different estimates of upwelling in ERA-Interim . The nudged simulations, especially the simulations RC1SD-base-09/10 without global mean temperature nudging, are closer to ERA-Interim data than the RC1-simulations without nudging. The strong influence of the SSTs in contrast to the relatively small influence of the vertical resolution on the tropical upward mass flux shows the strong SST effect on the stratospheric mass transport, known from earlier studies e.g.. The simulation RC2-oce-01 with interactively coupled ocean (Sect. ) shows a smaller mass flux than the other simulations and therefore fits best to ERA-Interim observations in the lower and middle stratosphere.

In the upper stratosphere the vertical resolution plays an increasing role. $F_{\text{trop}}$ is smaller in the simulations with 90 levels, which is in better agreement with ERA-Interim data. As the nudging does not extend up to 1 $\mathrm{hPa}$ (Sect. ), the RC1SD-simulations do not perform better than the other simulations.

Annual mean tropical upward mass flux (${10}^{\mathrm{8}}\mathrm{kg}{\mathrm{s}}^{-\mathrm{1}}$) at $70/74$ $\mathrm{hPa}$ (top), $\mathrm{9}/10$ $\mathrm{hPa}$ (middle), and 1 $\mathrm{hPa}$ (bottom). The simulation results (colours) are compared to ERA-Interim data (black).

[Figure omitted. See PDF]

Mean age of stratospheric air

Figure  shows the mean age of stratospheric air from the simulations (colours) and MIPAS data (black) averaged from 2002 to 2011 at around 50 $\mathrm{hPa}$ (20 $\mathrm{km}$ for MIPAS). The simulation results group according to their vertical resolution: the simulations with L47 show a lower mean AoA at all latitudes compared to the simulations with 90 levels, i.e. a stronger stratospheric transport. Results of the vertically higher-resolved set-ups are in better agreement with the MIPAS observations, the L90 results fit well with a mean AoA of around 3 $\mathrm{years}$ at mid-latitudes from MIPAS. Moreover, the nudged simulations perform slightly better compared to the MIPAS observations than the free-running simulations. In the tropics and at higher latitudes mean AoA from MIPAS is higher by about one-half and up to 2 years, respectively. However, note that AoA derived from ${\mathrm{SF}}_{\mathrm{6}}$ measurements by MIPAS is affected by the mesospheric sink of ${\mathrm{SF}}_{\mathrm{6}}$, which is not represented in the model. In particular in high latitudes, this effect can lead to an overestimation of AoA as derived from ${\mathrm{SF}}_{\mathrm{6}}$ of up to 2 $\mathrm{years}$ . A lower mean AoA and therefore a faster mass transport from model simulations compared to observations is well known from other CCMs . The simulation with coupled ocean lies in between the vertically lower and higher-resolved simulations and thereby fits better with observations from MIPAS than the L47 simulations with prescribed – observed or simulated – SSTs.

The time evolution of mean AoA from the model simulations at different stratospheric heights and latitudinal bands is depicted in Fig. . As MIPAS data are available for 10 years only, it is included as a mean over the years 2002 to 2011 with the corresponding standard deviation. Consistent with the rising trend in $F_{\text{trop}}$ (Sect. ), all simulations show a decrease in the mean AoA in the entire stratosphere. The implied BDC increase is the strongest in the lower stratosphere at mid- and high latitudes.

Mean age of stratospheric air ($\mathrm{years}$) at 47 $\mathrm{hPa}$$/$48 $\mathrm{hPa}$$/$20 $\mathrm{km}$, 90${}^{\circ }$ N–90${}^{\circ }$ S, averaged for the years 2002 to 2011. The simulation results (colours) are compared to MIPAS data (black).

[Figure omitted. See PDF]

The simulations with 90 levels show higher mean AoA at all stratospheric layers and latitudinal bands and therefore fit better with MIPAS observations. In the lower stratosphere (Fig. , top), the nudged simulations RC1SD-base-07 and RC1SD-base-10 are in best agreement with the observations, especially at mid-latitudes. The simulation RC2-oce-01 lies in between the L90 simulations and the vertically lower-resolved simulations. In contrast to the findings for the RC (Sect. ), the vertical resolution is much more important for the mean AoA than different SST data sets.

In the middle and upper stratosphere (Fig. , middle and bottom) the differences between the simulations become smaller and the mean AoA from observations gets considerably larger than from the simulations. One reason is the different derivation of the mean AoA: for EMAC an inert synthetic tracer is used, compared to the ${\mathrm{SF}}_{\mathrm{6}}$ tracer for the observations. ${\mathrm{SF}}_{\mathrm{6}}$ is inert in most parts of the stratosphere, but photochemically destroyed in the mesosphere. The decent of ${\mathrm{SF}}_{\mathrm{6}}$-poor air from the mesosphere into the stratosphere raises the calculated mean AoA . This might explain some of the deviations between EMAC and MIPAS data in the upper stratosphere.

Time series of mean age of stratospheric air ($\mathrm{years}$) at 74 $\mathrm{hPa}$$/$18 $\mathrm{km}$ (top), 9 $\mathrm{hPa}$$/$10 $\mathrm{hPa}$$/$30 $\mathrm{km}$ (middle), and 1 $\mathrm{hPa}$$/$50 $\mathrm{km}$ (bottom), averaged over the tropics (30${}^{\circ }$ S–30${}^{\circ }$ N, left column), mid-latitudes (30–60${}^{\circ }$ N/S, middle column), and high latitudes (60–90${}^{\circ }$ N/S, right column) in the annual mean. The simulation results (colours) are compared to MIPAS data (black), which is averaged for the years 2002 to 2011 with the corresponding standard deviation.

[Figure omitted. See PDF]

In summary, EMAC is able to simulate the Brewer–Dobson circulation in terms of the tropical upward mass flux and the mean age of stratospheric air in reasonable agreement with ERA-Interim reanalysis data and MIPAS observations, respectively. The simulations with 90 levels show a slower transport in the lower stratosphere than the simulations with the lower vertical resolution and are therefore in better agreement with observations. The simulation RC2-oce-01 with coupled ocean performs best concerning the lower and middle stratospheric RC. Mean AoA from this simulation lies in between the resolutions L47 and L90.

Tropospheric and stratospheric ozone

In this section we evaluate the ozone distributions of the simulations by comparing with observations. Data used as references are BSTCO Bodeker Scientific combined total column ozone database; for total column ozone (TCO), AURA Microwave Limb Sounder/Ozone Monitoring Instrument (MLS/OMI; ) for tropospheric and stratospheric partial column ozone (TPCO, SPCO), and the ozone sonde data set described by for ozone profiles. Since the data sets and simulations cover different time periods, the analyses are performed for the periods 1980–2011 and 2005–2011 for TCO and the partial columns, respectively, except for the RC1-aero-07 simulation which starts in 1990. The comparison with the ozone sonde data is based on the period 1995 to 2011. Time series are shown for all years that are available. Simulation data and observations are regridded to the coarsest common grid and represented on the same time axis for each comparison. To test the statistical significance we applied the paired $t$ test, based on observation–simulation pairs of values from each step of the respective time axes (diagnostic-specific). The portrait diagrams for overall mean bias metrics were adapted from the corresponding analyses of , and also share some routines with the ESMValTool .

Figure  provides an overview of the TCO bias in the different simulations with respect to the BSTCO database, globally and for different latitude bands. Bias and $t$ test calculations are based on annual data, spatially averaged for the corresponding latitude bands. Weighting considers grid cell area and the number of days per month. The null hypothesis of a zero bias can be rejected with a significance level of at least 90 $\mathrm{\%}$ for each tile (not shown). TCO is overestimated in all simulations, and the bias generally increases from north to south. Nudging reduces the bias, particularly in the SH high latitudes. We also note that the simulations with global mean temperature nudging (RC1SD-base-07, RC1SD-base-08) agree better with BSTCO reference data than the corresponding simulations without nudging wave zero of the temperature field (RC1SD-base-10, RC1SD-base-09).

Mean total column ozone bias (in %) between the simulations and the reference data set BSTCO for different latitude bands. The analyses are based on annual means of the years 1980 to 2011. The only exception is RC1-aero-07, which starts 1990.

[Figure omitted. See PDF]

Left column (a–e): climatological annual cycle of total column ozone (1980–2011, in $\mathrm{DU}$) from BSTCO observations (a) and selected simulations (RC1-base-07: b; RC2-base-04: c; RC1SD-base-07: d; RC2-oce-01: e). Middle column (f–i): difference between the simulations and the BSTCO observations (in $\mathrm{DU}$). Right column (j–l): difference (in $\mathrm{DU}$) between the corresponding simulations with high (b–e) and low vertical resolution (RC1-base-08: j; RC2-base-05: k; RC1SD-base-08: l). Statistical significance is tested with the paired $t$ test. Statistically significant changes on the 95 $\mathrm{\%}$ confidence level are coloured.

[Figure omitted. See PDF]

Figure  shows the climatological annual cycle of TCO for the BSTCO observations (Fig. a) and the high-resolution (L90MA) RC1 (Fig. b), RC2 (Fig. c) and RC1SD (Fig. d) simulations as well as the simulation with coupled ocean (RC2-oce-01, Fig. e). In general the seasonal cycle and the spatial distribution are well reproduced in the simulations with low TCO in the tropics and maximum values in the NH high latitudes in winter and spring. The minimum TCO values occur in SH spring in the SH polar region, which is referred to as the ozone hole. However, the temporal evolution of the ozone hole is not fully captured and the minimum TCO values are mostly underestimated, i.e. too much ozone is simulated. The differences (Fig. f–i) show that in SH spring the positive bias related to the underestimation of the ozone hole is the largest in the free-running RC1-base-07 (Fig. f) and RC2-oce-01 (Fig. i) simulations, whereas the smallest bias is found in the nudged simulations (Fig. h). Further analysis showed that the underestimation of the ozone hole in the free-running model simulation results from a too weak polar vortex, also apparent as too high temperatures in SH high latitudes in Fig. . We found that (planetary scale) wave fluxes in SH mid-latitudes are overestimated, and thus the polar vortex is too disturbed. The overestimation of SH spring TCO in the free-running simulations is smallest in the RC2 simulations with prescribed (simulated) SSTs/SICs (Fig. g). This is caused by a reduced (wave number 1) heat flux in the SH mid-latitudes for simulated compared to observed SSTs, which leads to a stronger polar vortex from August to October and lower temperatures inside the vortex (not shown). This favours the formation of polar stratospheric clouds and hence the heterogeneous ozone destruction.

Time series of annual mean total column ozone (in $\mathrm{DU}$) averaged over different latitude bands. BSTCO observations are shown in black. Note the different scales for the different regions.

[Figure omitted. See PDF]

The best agreement is found in the NH polar region in winter and spring in the free-running L90MA simulations (Fig. f, g). Comparing the L90MA and L47MA simulations (Fig. j–l), we find higher TCO values (and thus a larger bias) in the extra-tropics in the L47MA free-running simulations, but smaller values (and thus a smaller bias) in the nudged simulations.

Left column: climatological annual mean of observed total column ozone (BSTCO, top) as well as stratospheric (middle) and tropospheric (bottom) partial column ozone (AURA MLS/OMI), all in $\mathrm{DU}$. The tropospheric and stratospheric partial columns are integrated from the surface to the (on-line diagnosed, see text) tropopause and above the tropopause, respectively. The analyses cover the years 2005 to 2011. Middle column: differences between the RC1-base-07 simulation and the observations. Right column: same as middle column, but for the RC1SD-base-07 simulation. Statistically significant changes on the 95 $\mathrm{\%}$ confidence level are coloured.

[Figure omitted. See PDF]

The evolution of TCO between 1960 and 2100 is analysed for different latitude bands (Fig. ). Related to the increase of ODSs during the 20th century and their regulation since the 1990s e.g., TCO shows a negative trend until the 2000s and a positive trend in the 21st century in the extra-tropical regions. The decrease in the nudged simulations is comparable to the observations, whereas it is slower in the free-running simulations. In the future, all simulations project an increase of extra-tropical TCO with a return to 1980 levels until 2100. This positive trend is consistent with earlier CCM studies e.g.. At SH mid- and high latitudes the return to 1980 levels is delayed in the L90MA simulations compared to the L47MA simulations, which is linked to the larger ozone loss in the 1990s (see Fig. S35). In the tropics, TCO slightly increases in the first half of the 21st century in all RC2 simulations, and decreases afterwards. This is related to a strengthening of the BDC (see Sect. ) during the 21st century, which leads to an enhanced export of lower stratospheric ozone and thus counteracts the chemically induced ozone increase. To better understand the differences between the simulations and the observations, the partial column ozone for the troposphere and the stratosphere is analysed (Fig. ). Regarding the free-running L90MA RC1 simulation, a significant positive TCO bias to the observations is found in the tropics and the SH. The main contribution to the overestimation of TCO in the tropics results from tropospheric ozone, while in the SH the stratospheric ozone bias is larger. Including global mean temperature in the nudging leads to a better agreement in the tropical tropospheric ozone column and the SH stratospheric ozone column, but enlarges the positive bias of stratospheric ozone at NH mid-latitudes. For all simulations the mean bias of tropospheric column ozone to the observation is shown in Fig. . It is positive in all regions and simulations see also. Taking into account the corrected road traffic emissions (RC1SD-base-10a), the tropospheric partial column ozone increases by up to 3 % in the tropics and thus also the bias is increased (not shown). Furthermore, the comparison of the two nudging methods shows again that the mean bias is mostly reduced, if nudging of the global mean temperature is included. For this, at least three effects are potentially responsible: (1) as the temperature mean nudging increases overall the tropospheric temperature by up to 4 $\mathrm{K}$ (see Sect. ), the temperature-dependent reaction kinetics is altered, pushing the chemistry into a new equilibrium state e.g.. (2) The TPCO is reduced due to a decreased tropopause height, if mean temperature is nudged. Indeed, the lower tropopause height (higher tropopause pressure by up to $+$10 $\mathrm{hPa}$) of RC1SD-base-07 (with $\stackrel{\mathrm{‾}}{T}$ nudging) compared to RC1SD-base-10 (without $\stackrel{\mathrm{‾}}{T}$ nudging) reduces the TPCO polewards of 40${}^{\circ }$ latitude by up to $-$0.4 $\mathrm{DU}$ (Fig. S38). Around 30${}^{\circ }$ N/S, the corresponding shift is up to $-$5 $\mathrm{hPa}$ increasing the TPCO by up to 1.5 $\mathrm{DU}$. The resulting global effect of this geometrical tropopause shift (analysed for the 11-year period 2000 to 2010) is $-$0.1 $\mathrm{DU}$. This effect is smaller, however, than the influence of the tropopause definition: using the on-line diagnosed tropopause (WMO definition equatorward of 30${}^{\circ }$ latitude; iso-surface of 3.5 $\mathrm{PVU}$ potential vorticity poleward of 30${}^{\circ }$ latitude) of the submodel TROPOP instead of the WMO tropopause partially compensates the high bias in EMAC (poleward of 30${}^{\circ }$ latitude: up to $-$4 $\mathrm{DU}$; global average (2000–2010): $-0.4$ $\mathrm{DU}$; see Fig. S39) when comparing to MLS/OMI, where TPCO is based on the WMO tropopause (but another temperature field). (3) The vertical temperature gradient is reduced by the global mean temperature nudging, thus reducing the convective activity and with it the production of $\mathrm{NO}$${}_x$ from lightning. Indeed, this is the case (see Sect.  and Fig. ): global mean temperature nudging reduces the lightning $\mathrm{NO}$${}_x$ production from 4.5 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$ (RC1SD-base-10) to 1.8 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$ (RC1SD-base-07), i.e. by 2.7 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$. Using the results of their Fig. 3b this explains a reduction of the ozone column by about 5.4 $\mathrm{DU}$ (i.e. approximately 2 $\mathrm{DU}\left(\mathrm{Tg}\right(\mathrm{N}){\mathrm{a}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$). This is based on their finding that the average ozone production efficiency by $\mathrm{NO}$${}_x$ from lightning is 5 times higher compared to other $\mathrm{NO}$${}_x$ sources (except for aviation), and that 70 $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$ on average account for about 26.5 $\mathrm{DU}$ ozone (i.e. 0.4 $\mathrm{DU}\left(\mathrm{Tg}\right(\mathrm{N}){\mathrm{a}}^{-\mathrm{1}}{)}^{-\mathrm{1}}$ on average for all sources). This dominating effect of the modified lightning $\mathrm{NO}$${}_x$ emissions is further corroborated by the changed ozone production and loss rates, which are altered by the global mean temperature nudging mostly in the tropics, where the lightning $\mathrm{NO}$${}_x$ production peaks (Fig. S40, Appendix ).

Same as Fig. , but for the tropospheric partial column ozone. The reference data set is AURA MLS/OMI. All tiles are based on annual mean values of the years 2005 to 2011 and statistically significant on the 90 % confidence level. Here, the WMO tropopause definition has been used to be as consistent as possible with the reference data set.

[Figure omitted. See PDF]

Additional effects, which are however more difficult to quantify, are direct effects on ozone by altered convection, by altered mixing, or by modified stratosphere–troposphere exchange.

In contrast to TCO, TPCO increases from 1960 to the second half of the 21st century and slightly decreases afterwards (Fig. S36). This development is linked to the prescribed emissions of tropospheric ozone precursors, in particular methane, according to the RCP6.0 scenario .

Same as Fig. , but for the ozone volume mixing ratio at 6 selected pressure levels (700, 500, 250, 100, 50, 10 $\mathrm{hPa}$). The ozone sonde data set is used as reference. All tiles are based on the annual cycle, calculated from monthly data of the years 1995 to 2011. Results of a paired $t$ test for the probability to reject the null hypothesis of equal means is shown on the right.

[Figure omitted. See PDF]

Vertically resolved ozone is compared to ozone sonde data . This analysis has been adapted from , also sharing some routines with the ESMValTool . Simulations were sampled at the ozone sonde locations, binned according to latitude ranges, and averaged with equal weights. All simulation output time steps (every 10 $\mathrm{h}$) from 1995 to 2011 contributed to the calculation of annual cycle data. Ozone sonde data are from the same period, but the annual cycle is based on less samples. Thus, data are co-located in space, but not necessarily in time. Most tiles are statistically significant (90 $\mathrm{\%}$ level), but each $t$ test is based on 12 data pairs only. Figure  shows the mean bias and corresponding significance for different pressure levels and latitude bands. The comparisons to ozone sonde profile data generally confirm the findings for total and partial column ozone discussed above: There is mostly a high bias, increasing from north to south, with nudged simulations performing best. The annual cycle of ozone volume mixing ratio is generally well reproduced in the stratosphere (except in the SH polar region) as well as in the troposphere at mid-latitudes (Fig. S37). Some differences occur, however, in the upper troposphere in the tropics and at high latitudes. At the 250 $\mathrm{hPa}$ level the bias is strongly negative in northern and southern high latitudes, reversing to positive towards the tropics. The contrast of the bias in 250 $\mathrm{hPa}$ and neighbouring levels is largest in the free-running simulations, indicating a dynamical problem in the tropopause region. We also note that the positive bias in the lower stratosphere (100 and 50 $\mathrm{hPa}$) is considerably larger in the SH, which is linked to the underestimated ozone hole.

Summary and conclusions

With the chemistry–climate model EMAC (version 2.51), we performed a set of reference simulations as recommended by the Chemistry–Climate Model Initiative (CCMI): hindcast simulations (1950–2011) without and with specified dynamics and combined hindcast and projection simulations (1950–2100) based on the RCP6.0 scenario. We performed simulations at T42 spectral resolution with two different vertical resolutions, L90MA and L47MA with 90 and 47 hybrid model layers between the surface and approximately 80 $\mathrm{km}$ altitude (in the MA), respectively. One simulation (1950–2100) in T42L47MA was performed with an interactively coupled ocean model and set up based on an extensive spin-up procedure. For the simulation with specified dynamics, two different Newtonian relaxation (nudging) set-ups, both using ERA-Interim reanalysis data, have been applied: either excluding or including the global mean temperature. Additional hindcast simulations have been performed with additional on-line calculated tropospheric aerosol, with and without coupling to the cloud processes. All simulations have been equipped with comprehensive on-line diagnostics.

The manuscript describes briefly the EMAC model updates and in detail the different model set-ups, including some unintended deviations from the CCMI recommendations and corresponding sensitivity simulations. The description also includes the applied on-line diagnostics and an analysis of the on-line calculated source (primary emissions) and sink (dry and wet deposition, aerosol sedimentation) terms, and is meant as data set description and reference for further analysis of the close to 2 $\mathrm{PetaByte}$ comprising data set.

First analyses presented here focus mainly on an intercomparison between the different simulations from a global perspective. The precipitation patterns of all simulations are within the range of results of comparable models, and the coupled atmosphere–ocean model shows the largest deviations from observations. The temperature distributions in the nudged mode (without global mean temperature nudging) show a tropospheric and stratospheric cold bias of up to 4 $\mathrm{K}$. In free-running mode with prescribed SSTs and independent of the vertical resolution, the cold bias is largest in the UTLS and the tropical stratosphere, whereas a pronounced warm bias appears in the SH extratropics-to-polar transition region (50–70${}^{\circ }$ S) above 140 $\mathrm{hPa}$ and in the SH polar region between 140 and 30 $\mathrm{hPa}$. This warm bias is significantly reduced in the simulation with coupled ocean model.

For the stratospheric mean age of air, the vertical resolution has the largest impact: the mean age of air is on average by about 1 year younger in L47MA compared to L90MA in the lower stratosphere at mid- and high latitudes. The nudged simulations with 90 layers show the best agreement with observations. The simulation RC2-oce-01 performs best concerning the lower and middle stratospheric residual circulation. For the tropical upward mass flux the impact of the SST is larger compared to that of the vertical resolution.

Independent of the model set-up and the model resolution, all free-running simulations show a positive bias of the total ozone column of up to 40 $\mathrm{DU}$ increasing from about 5 $\mathrm{DU}$ over the north pole to the maximum values in the SH polar region. Nudging reduces this bias, in particular in the SH polar region. Independent of the set-up, the ozone bias in the tropics is mainly caused by overestimated tropospheric ozone, whereas in the extratropics and polar regions the contribution of stratospheric ozone increases with latitude.

Including the global mean temperature in the nudging procedure increases the tropospheric temperature by up to 4 $\mathrm{K}$ throughout the year and weakens the vertical gradient between the surface and the tropopause (which is also shifted). This alters the convective activity with a large impact on the production of $\mathrm{NO}$${}_x$ by lightning, at least in the applied updraft velocity-based scheme. With this, we confirm an earlier rule of thumb estimate deriving an increase of about 2 $\mathrm{DU}$ in tropospheric total ozone column per $\mathrm{Tg}\left(\mathrm{N}\right){\mathrm{a}}^{-\mathrm{1}}$ lightning $\mathrm{NO}$${}_x$ production. The reduced lightning $\mathrm{NO}$${}_x$ production with global mean temperature nudging is partly compensated by an increased $\mathrm{NO}$${}_x$ release from soil, due to increased soil temperatures. Whereas $\mathrm{NO}$${}_x$ emissions from soil increase with the global mean temperature nudging, the effect on isoprene emissions from the biosphere is less pronounced. Emissions from the ocean ($\mathrm{DMS}$, ${\mathrm{C}}_{\mathrm{5}}{\mathrm{H}}_{\mathrm{8}}$, $\mathrm{Br}$), in contrast, are reduced by nudging, and even more reduced, if global mean temperature is nudged. The oceanic uptake of ${\mathrm{CH}}_{\mathrm{3}}\mathrm{OH}$ is likewise reduced with (global mean temperature) nudging. Deposition fluxes through scavenging and dry deposition of ozone, nitrate, $\mathrm{NO}$${}_x$, and sulfite are reduced by the global mean temperature nudging at both vertical resolutions, other deposition fluxes are less affected. At least for ozone, nitrate, and $\mathrm{NO}$${}_x$ this is consistent to the likewise-reduced lightning $\mathrm{NO}$${}_x$ production and the reduced ${\mathrm{O}}_{\mathrm{3}}$ burden.

The increased tropospheric temperature through the global mean temperature nudging slightly reduces the tropospheric methane lifetime towards $\mathrm{OH}$, this effect is larger in simulations with the coarser vertical resolution (L47MA). Overall, i.e. averaged over all simulations, the simulated $\mathrm{OH}$ lifetime of methane is $8.0\pm 0.6$ $\mathrm{years}$ in the period 2000–2004. This is potentially too short, indicated by the consistently underestimated ${\mathrm{CH}}_{\mathrm{4}}$ mixing ratios in the NH UTLS, i.e. compared to observations at CARIBIC flight levels. Note that ${\mathrm{CH}}_{\mathrm{4}}$ has been prescribed (based on observations) by Newtonian relaxation at the lower boundary.

The simulation results will be made publicly available (see next section) for further in-depth analyses. For intercomparison with observations, we recommend to use the results of the nudged simulation with all corrections, i.e. RC1SD-base-10a. The results of RC1SD-base-07 and RC1SD-base-08 should be used with caution, due to the large impact of the global mean temperature nudging, for which no specific parameter re-optimization for the radiation balance has been undertaken yet. Such an optimization will certainly alter the hydrological cycle, i.e. clouds and convection, and with it also the lightning $\mathrm{NO}$${}_x$ production. Studies for which the specified dynamics (nudging) is not desired, e.g. on trends and frequency distributions, are best based on the results of the free-running simulations with 90-level discretization. Nevertheless, any intercomparison to those with 47 levels is also desirable, in particular since the simulation with coupled ocean model was performed with 47 levels in the atmosphere. Last, but not least, for further analyses on aerosol and aerosol–cloud effects, only RC1-aero-07 (from 1991 onwards) and RC1-aecl-02 (from 1966 onwards) should be used, respectively.

Code and data availability

The Modular Earth Submodel System (MESSy) is continuously further developed and applied by a consortium of institutions. The usage of MESSy and access to the source code is licensed to all affiliates of institutions, which are members of the MESSy Consortium. Institutions can become a member of the MESSy Consortium by signing the MESSy Memorandum of Understanding. More information can be found on the MESSy Consortium Web-site (http://www.messy-interface.org).

The data of the simulations described above will be made available in the Climate and Environmental Retrieval and Archive (CERA) database at the German Climate Computing Centre (DKRZ; http://cera-www.dkrz.de/WDCC/ui/Index.jsp). The corresponding digital object identifiers (doi) will be published on the MESSy consortium web-page (http://www.messy-interface.org). A subset of the data of those simulations covering consistently the requested time periods (1960–2010 for RC1, and 1960–2099 for RC2) will be submitted to the BADC database for the CCMI project.

Additional on-line diagnostics

TBUDGET

This submodel analyses the contribution of different production cycles (e.g. associated with different source gases) to selected tracers type-2 family, see. This is done by defining additional diagnostic tracers to store corresponding production and loss rates. These are specified by namelist entries together with the total tracer and the loss of the total tracer. For a tracer $C$, a component-tracer $C_i$ is calculated by $${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(C_i\right)=P_i+L\times \left({\displaystyle \frac{C_i}{C}}\right),$$ where $L$ is the destruction rate of $C$ and $P_i$ the production of tracer $C_i$. A numerical correction (scaling) of $C_i$ to ensure that the sum (over $i$) of all $C_i$ equals the total tracer (family) $C$ is applied automatically. In the simulations this has been applied to distinguish short- from long-lived halogenated species to assess the influence of VSLS (very short-lived species) on the ozone budget; see Figs. S40–S43: Type-2 tracer families $\mathrm{Br}$${}_y$ and $\mathrm{Cl}$${}_y$ are defined in the CTRL_FAMILY namelist of submodel TRACER, submodel SCALC is used to sum the loss rates $\mathrm{LossBr}$ and $\mathrm{LossCl}$ from their individual tendencies calculated by SCAV, and in the CPL namelist of TBUDGET the corresponding diagnostic tracers $\mathrm{BrS}$, $\mathrm{BrL}$, $\mathrm{ClS}$, and $\mathrm{ClL}$ for short (S)- and long (L)-lived halogen compounds, respectively, are defined for integration according to Eq. (). In addition, tracers for reactive halogen production from L- and S-lived halogenated compounds, respectively, have been defined within the MECCA chemical mechanism (see ESCiMo_MECCA_mechanism.pdf in the Supplement): $\mathrm{ProdLBr}$, $\mathrm{ProdSBr}$, $\mathrm{ProdLCl}$, $\mathrm{ProdLBr}$. Note that these diagnostic tracers contain the accumulated (over time) rates.

O3ORIG

The ozone origin diagnostics of has been applied to define 14 diagnostic ozone tracers ($\mathrm{O}\mathrm{3}\mathrm{ONHTS}$, $\mathrm{O}\mathrm{3}\mathrm{OTRTS}$, $\mathrm{O}\mathrm{3}\mathrm{OSHTS}$, $\mathrm{O}\mathrm{3}\mathrm{ONPLS}$, $\mathrm{O}\mathrm{3}\mathrm{ONMLS}$, $\mathrm{O}\mathrm{3}\mathrm{OTRLS}$, $\mathrm{O}\mathrm{3}\mathrm{OTRMS}$, $\mathrm{O}\mathrm{3}\mathrm{OSMLS}$, $\mathrm{O}\mathrm{3}\mathrm{OSPLS}$, $\mathrm{O}\mathrm{3}\mathrm{ONPUS}$, $\mathrm{O}\mathrm{3}\mathrm{ONMUS}$, $\mathrm{O}\mathrm{3}\mathrm{OTRUS}$, $\mathrm{O}\mathrm{3}\mathrm{OSMUS}$, $\mathrm{O}\mathrm{3}\mathrm{OSPUS}$). Each ozone tracer field is produced only in one specific region: NHTS (Northern Hemisphere troposphere), TRTS (tropical troposphere), SHTS (Southern Hemisphere troposphere), NPLS (northern polar lower stratosphere), NMLS (northern mid-latitudes lower stratosphere), TRLS (tropical lower stratosphere), TRMS (tropical middle stratosphere), SMLS (southern mid-latitudes lower stratosphere), SPLS (southern polar lower stratosphere), NPUS (northern polar upper stratosphere), NMUS (northern mid-latitudes upper stratosphere), TRUS (tropical upper stratosphere), SMUS (southern mid-latitudes upper stratosphere), and SPUS (southern polar upper stratosphere). The sum of all 14 ozone origin tracers is the full ozone field. The definition of the regions in terms of height and latitude is illustrated in Fig. S45.

CONTRAIL

The submodel CONTRAIL diagnoses the potential contrail coverage (variable potcov), describing the fractional area in which contrails can form and persist according to the Schmidt–Appleman criterion . In addition, the potential contrail cirrus coverage (variable b_cc) is diagnosed, taking into account regions where contrails can persist once they have been formed. Both variables were calculated according to and were output (channel “contrail_gp”) as 10-hourly global snapshots.

Specific sampling of model data

With the submodel SCOUT Selectable Column OUTput;Sect. 5.2 hourly output of vertical profiles (i.e. the model column) has been sampled during all simulations at 49 stations of the NOAA/ESRL network, 16 stations of the SHADOZ network, and 11 stations of the WOUDC. The complete namelist for reference is shown in Fig. S47.

With the submodel SORBIT Sampling along ORBITs,Sect. 5.4 data along sun-synchronous orbits of eight different satellites have been output during all simulations. The complete namelist for reference is shown in Fig. S48.

With the submodel VISO iso-surfaces and maps;Sect. 5.1 isentropes of 340, 380, and 420 $\mathrm{K}$ and iso-surfaces of 2, 3, and 4 $\mathrm{PVU}$ potential vorticity (PV) are defined. On the PV iso-surfaces, pressure, potential temperature, and temperature are mapped; pressure and PV are mapped on the isentropes. Further temperature, geopotential, and potential temperature are mapped on the tropopause (TP), and temperature and pressure on the planetary boundary layer height (PBLH). See Fig. S46 for the complete namelist. PBLH and TP are calculated by the submodel TROPOP, TP is defined according to the WMO equatorward of 30${}^{\circ }$ latitude, and as 3.5 $\mathrm{PVU}$ iso-surface poleward. All corresponding surfaces and maps are output as 10-hourly snapshots.

With the submodel S4D sampling in 4 dimensions;Sect. 5.3 data along (available) tracks of several research platforms have been sampled during the simulations with specified dynamics (RC1SD-base). The complete namelist for reference is shown in Fig. S48.

In addition to the 10-hourly global model output, 6-hourly global snapshots have been written (channel “6h”) for the horizontal wind velocity components (um1, vm1), temperature (tm1), specific humidity (qm1), the vertical velocity (etadot, vervel), and the geopotential (geopot).

The submodel SATSIMS was used for the on-line calculation of cloud optical properties comparable to the International Satellite Cloud Climatology Project (ISCCP) data set with the ISCCP cloud simulator . This uses cloud optical thickness and cloud top pressure to derive cloud types according to the ISCCP classification.

Monthly mean values (on-line calculated, channel “mm”) and corresponding standard deviations have been output (in addition) for the tracers ${\mathrm{CO}}_{\mathrm{2}}$, ${\mathrm{O}}_{\mathrm{3}}$, ${\mathrm{CH}}_{\mathrm{4}}$, ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$, ${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$, ${\mathrm{CF}}_{\mathrm{2}}{\mathrm{Cl}}_{\mathrm{2}}$, ${\mathrm{CFCl}}_{\mathrm{3}}$, ${\mathrm{HNO}}_{\mathrm{3}}$, ${\mathrm{HNO}}_{\mathrm{3}}\mathrm{\_}\mathrm{nat}*$, $\mathrm{OH}$, ${\mathrm{O}}^{\mathrm{1}}\mathrm{D}$, $\mathrm{Cl}$, and the photolysis rate ${\mathrm{J}}_{\mathrm{CH}\mathrm{4}}$. Note that this additional output can directly be used to set-up subsequent EMAC simulations in quasi-chemistry-transport model (QCTM) mode .

Additional diagnostic tracers

The submodel DRADON diagnostic Radon,Sect. 6.1 has been used to simulate ${}^{222}\mathrm{Rn}$ and ${}^{210}\mathrm{Pb}$ with a constant ${}^{222}\mathrm{Rn}$ source of 10 000 $\mathrm{atoms}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ over ice-free land (zero elsewhere). Aerosol properties ($\mathit{\sigma }$ $=$ 2, radius $=$ $\mathrm{5}\times {10}^{-\mathrm{7}}$ $\mathrm{m}$) have been assigned to ${}^{210}\mathrm{Pb}$ by the submodel PTRAC (via a dummy $\mathrm{AEROSOL}$ tracer), so that the loss processes (submodels SCAV, SEDI and DDEP) are simulated realistically.

For additional analyses of the ozone budget, accumulated ozone production and loss rates (diagnostic tracers) have been defined within the MECCA chemical mechanism ($\mathrm{ProdO}\mathrm{3}$, $\mathrm{LossO}\mathrm{3}$, redirected for output to channel “o3orig”), as well as contributions by individual reactants or reactant families: $\mathrm{LossO}\mathrm{3}\mathrm{Br}$, $\mathrm{LossO}\mathrm{3}\mathrm{Cl}$, $\mathrm{LossO}\mathrm{3}\mathrm{H}$, $\mathrm{LossO}\mathrm{3}\mathrm{N}$, $\mathrm{LossO}\mathrm{3}\mathrm{N}\mathrm{2}$, $\mathrm{LossO}\mathrm{3}\mathrm{O}$, $\mathrm{LossO}\mathrm{3}\mathrm{R}$, $\mathrm{LossO}\mathrm{3}\mathrm{Hn}$, $\mathrm{LossO}\mathrm{3}\mathrm{Nn}$, $\mathrm{LossO}\mathrm{3}\mathrm{Cln}$ (channel “tr_o3_bud”); $\mathrm{LossHO}\mathrm{2}$, $\mathrm{LossO}\mathrm{1}\mathrm{D}$, $\mathrm{LossOH}$ (channel “tr_o3_tbud”). All were output as 10-hourly global snapshots. Note that these diagnostic tracers contain the accumulated (over time) rates.

An additional set of tracers, $({\mathrm{CFCl}}_{\mathrm{3}}{)}_{\mathrm{c}}$, $({\mathrm{CF}}_{\mathrm{2}}{\mathrm{Cl}}_{\mathrm{2}}{)}_{\mathrm{c}}$, $({\mathrm{CH}}_{\mathrm{3}}{\mathrm{CCl}}_{\mathrm{3}}{)}_{\mathrm{c}}$, $({\mathrm{N}}_{\mathrm{2}}\mathrm{O}{)}_{\mathrm{c}}$, $({\mathrm{CF}}_{\mathrm{2}}\mathrm{ClBr}{)}_{\mathrm{c}}$, and $({\mathrm{CF}}_{\mathrm{3}}\mathrm{Br}{)}_{\mathrm{c}}$, have been included in the MECCA chemical mechanism: they were relaxed (with submodel TNUDGE, relaxation time constant $\mathit{\tau }=\mathrm{3}$ $\mathrm{h}$) towards a constant mixing ratio of 100 $\mathrm{pmol}{\mathrm{mol}}^{-\mathrm{1}}$ at the surface and reacted in the same way as the corresponding species without index c (see ESCiMo_MECCA_mechanism.pdf in the Supplement).

As proposed by CCMI , several additional diagnostic tracers were included as listed in Table .

Additional diagnostic tracers. Acronyms in parentheses denote the submodels simulating the indicated processes. Corresponding data files are listed in Table E1.

Tracer	Defining	Description
	submodel
$\mathrm{SF}\mathrm{6}$	PTRAC${}^{\text{a}}$	Newtonian relaxation (TNUDGE)${}^{\text{c}}$ towards time series based on observations (see Fig. E6)
$\mathrm{AOA}$	PTRAC	Age of air tracer; Newtonian relaxation (TNUDGE)${}^{\text{c}}$ at lowest model layer towards linearly in time increasing mixing ratio
$\mathrm{SF}\mathrm{6}\mathrm{\_}\mathrm{AOA}$	PTRAC	Age of air tracer; Newtonian relaxation (TNUDGE)${}^{\text{c}}$ at lowest model layer towards a latitude dependent, linearly in time increasing mixing ratio
$\mathrm{SF}\mathrm{6}\mathrm{\_}\mathrm{AOAc}$	PTRAC	Age of air tracer; Newtonian relaxation (TNUDGE)${}^{\text{c}}$ at lowest model layer towards a linearly in time increasing mixing ratio
$\mathrm{SF}\mathrm{6}\mathrm{\_}\mathrm{CCMI}$${}^{\text{e}}$	PTRAC	Emissions (OFFEMIS) according to EDGAR v4.2 database (see Fig. E55)
${\mathrm{SO}}_{\mathrm{2}}\mathrm{t}$${}^{\text{f,g}}$	PTRAC	anthropogenic emissions (OFFEMIS) as ${\mathrm{SO}}_{\mathrm{2}}$, wet removal (SCAV) as ${\mathrm{SO}}_{\mathrm{2}}$
$\mathrm{NH}\mathrm{\_}05$${}^{\text{d}}$	TREXP${}^{\text{b}}$	Newtonian relaxation (TNUDGE)${}^{\text{j}}$ towards surface layer mixing ratio (100 $\mathrm{nmol}{\mathrm{mol}}^{-\mathrm{1}}$) over 30–50${}^{\circ }$ N, 5-day exponential decay ($e$-folding time $\mathit{\tau }=4.32\times {10}^{-\mathrm{5}}$ $\mathrm{s}$)
$\mathrm{NH}\mathrm{\_}50$${}^{\text{d}}$	TREXP	As $\mathrm{NH}\mathrm{\_}05$, but 50-day exponential decay
$\mathrm{NH}50\mathrm{W}$${}^{\text{d,g}}$	TREXP	As $\mathrm{NH}\mathrm{\_}50$, but additional wet removal (SCAV) as ${\mathrm{HNO}}_{\mathrm{3}}$
$\mathrm{AOA}\mathrm{\_}\mathrm{NH}$${}^{\text{d}}$	TREXP	Newtonian relaxation (TNUDGE)${}^{\text{j}}$ towards zero mixing ratio at surface level between 30 and 50${}^{\circ }$ N; uniform fixed source (OFFEMIS) of $3.171\times {10}^{-08}$ ($=$ years s${}^{-\mathrm{1}}$) $\mathrm{molecules}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$
$\mathrm{ST}80\mathrm{\_}25$${}^{\text{d}}$	TREXP	Newtonian relaxation (TNUDGE)${}^{\text{j}}$ towards 200 $\mathrm{nmol}{\mathrm{mol}}^{-\mathrm{1}}$ above ca. 80 $\mathrm{hPa}$ (level 61 in L90MA, level 20 in L47MA); 25-day exponential decay
$\mathrm{CO}\mathrm{\_}25$${}^{\text{h}}$	TREXP	Emitted (OFFEMIS) as anthropogenic $\mathrm{CO}$; 25-day exponential decay
$\mathrm{CO}\mathrm{\_}50$${}^{\text{h}}$	TREXP	Emitted (OFFEMIS) as anthropogenic $\mathrm{CO}$; 50-day exponential decay
$\mathrm{AOA}\mathrm{\_}\mathrm{CCMI}$${}^{\text{d}}$	TREXP	Newtonian relaxation (TNUDGE)${}^{\text{j}}$ towards zero mixing ratio at surface layer; uniform fixed source (OFFEMIS) as $\mathrm{AOA}\mathrm{\_}\mathrm{NH}$
${\mathrm{O}}_{\mathrm{3}}\left(\mathrm{s}\right)$${}^{\text{d,i}}$	MECCA	Stratospheric ozone tracer; Newtonian relaxation (TNUDGE)${}^{\text{j}}$ towards ${\mathrm{O}}_{\mathrm{3}}$ in the stratosphere; destroyed in the troposphere as ${\mathrm{O}}_{\mathrm{3}}$; the corresponding loss rate tracer $\mathrm{LO}\mathrm{3}\left(\mathrm{s}\right)$ is a qualitative measure for the troposphere to stratosphere exchange of ozone

${}^{\text{a}}$ ; ${}^{\text{b}}$  (2010; Tracer Release EXPeriment; Sect. 6.3); ${}^{\text{c}}$ relaxation time constant $\mathit{\tau }=\mathrm{3}$ $\mathrm{h}$; ${}^{\text{d}}$ according to CCMI (; Sect. 4.2); ${}^{\text{e}}$ dissenting from ${}^{\text{d}}$ without interpolation of emission time series to monthly values; ${}^{\text{f}}$ dissenting from ${}^{\text{d}}$ with transient anthropogenic emissions; ${}^{\text{g}}$ see ESCiMo_SCAV_mechanism.pdf in the Supplement; ${}^{\text{h}}$ dissenting from ${}^{\text{d}}$ with seasonal cycle of emission; ${}^{\text{i}}$ see ESCiMo_MECCA_mechanism.pdf in the Supplement; ${}^{\text{j}}$ relaxation time constant is equal to model time step length.

The Supplement related to this article is available online at doi:10.5194/gmd-9-1153-2016-supplement.

Acknowledgements

M. Righi and J. Hendricks (Deutsches Zentrum für Luft- und Raumfahrt (DLR), Institut für Physik der Atmosphäre, Oberpfaffenhofen, Germany) are acknowledged for providing the tropospheric aerosol climatology used for the base simulations.

We thank G. Bodeker of Bodeker Scientific for providing the combined total column ozone database.

We thank G. Stiller and F. Haenel (both at Karlsruhe Institute of Technology) for providing the MIPAS data and the support.

We thank all members of the IAGOS – CARIBIC team, who are not co-authors, especially P. van Velthoven (KNMI, the Netherlands), who provided the meteorological support. The collaboration with Lufthansa and Lufthansa Technik and the financial support from the German Ministry for Education and Science (grant 01LK1223C) are gratefully acknowledged. The data analysed in this paper can be accessed by signing the CARIBIC data protocol to be downloaded at http://www.caribic-atmospheric.com/.

We acknowledge the use of GPCP precipitation data products provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their website at http://www.esrl.noaa.gov/psd/ and the use of ERA-Interim precipitation data from ECMWF (synoptic monthly means of total precipitation; http://apps.ecmwf.int/datasets/data/interim-full-mnth/).

We acknowledge the World Climate Research Programme's Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the HadGEM climate modelling group for producing and making available their model output. For CMIP the US Department of Energy's Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led development of software infrastructure in partnership with the Global Organization for Earth System Science Portals.

We thank A. Voigt, S. Rast, and U. Schulzweida (all Max-Planck-Institute for Meteorology, Hamburg) for providing the modified ECHAM5 nudging code using netCDF input format and the corresponding preprocessing tools based on cdo (https://code.zmaw.de/projects/cdo).

The model simulations have been performed at the German Climate Computing Centre (DKRZ) through support from the Bundesministerium für Bildung und Forschung (BMBF). DKRZ and its scientific steering committee are gratefully acknowledged for providing the HPC and data archiving resources for this consortial project ESCiMo (Earth System Chemistry integrated Modelling). We also thank K. Fieg from DKRZ for servicing us for this project.

Parts of this work were produced by the DFG Research Unit FOR 1095 “Stratospheric Change and its Role for Climate Prediction” (SHARP) under the grant numbers LA 1025/13-2, LA 1025/14-2, LA 1025/15-2, JO 757/1-2, DA 233/3-1, DA 233/3-2, DA 233/4-1, DA 233/4-2, the BMBF project “Role Of the MIddle atmosphere in Climate: Quantification of Uncertainties of Solar-Induced Climate variability” (ROMIC-SOLIC, grant number 01LG1219B), the BMBF project “Mittelfristige Klimaprognose, The role of the stratosphere for decadal climate prediction” (MiKlip STRATO, grant number FKZ 01LP1168B), and the EU-project (FP7-ENV-2013) “Stratospheric and upper tropospheric processes for better climate predictions” (StratoClim, grant number: 603557).

We used the NCAR Command Language (NCL) for data analysis and to create some of the figures of this study. NCL is developed by UCAR/NCAR/CISL/TDD and available on-line: http://dx.doi.org/10.5065/D6WD3XH5. The authors wish further to acknowledge use of the Ferret program for analysis and graphics in this paper. Ferret is a product of NOAA's Pacific Marine Environmental Laboratory. Information is available at http://ferret.pmel.noaa.gov/Ferret/. Last but not least, we thank Veronika Eyring (DLR) for valuable comments on the first draft, as well as Fiona O'Connor (editor) and two anonymous referees for their constructive reviews. The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association. Edited by: F. O'Connor