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

. 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 speciﬁed 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 quantiﬁed and a ﬁrst 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 green-P.Jöckel et al.: ESCiMo with MESSy v. 2.51 house 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 highquality, 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 (Jöckel et al., 2006).For the first time, some of the Earth system models (ESMs) with interactive oceans participating in the Coupled Model Intercomparison Project Phase 5 (CMIP5; Taylor et al., 2011) had interactive chemistry (Eyring et al., 2013a;Flato et al., 2013).The World Meteorological Organization -United Nations Environment Programme (WMO/UNEP) "Scientific Assessment of Ozone Depletion: 2010" (World Meteorological Organisation, 2011) 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 (SPARC, 2010), that stratosphere-resolving CCMs should continue to evolve towards more comprehensive, self-consistent stratospheretroposphere CCMs.These developments provide a pathway for including better representation of stratospheretroposphere, chemistry-climate, and atmosphere-ocean coupling in CCMs and ESMs used for more robust predictions of future stratospheric ozone, climate change, and mutual influences (Eyring et al., 2013b).
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. 2 documents briefly (mainly for the users of it) the updates of the MESSy and the EMAC model since Jöckel et al. (2010), Sect. 3 describes the model set-ups, and Sect. 4 highlights some first analyses.Section 5 provides a summary and some first conclusions, followed by a section providing information about the code and data availability.Last but not least, Appendix A 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 multiinstitutional computer codes to the core atmospheric model, i.e. the 5th generation European Centre Hamburg general circulation model (ECHAM5; Roeckner et al., 2006).Since the publication of Jöckel et al. (2010), 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 (Kerkweg and Jöckel, 2015).The rediscretization via NCREGRID, as described by Jöckel (2006), is now part of those.Based on this advanced infrastructure, the submodels OFFLEM, ONLEM (both Kerkweg et al., 2006b), and DRYDEP (Kerkweg et al., 2006a) 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 (R. Sander et al., 2011) and further revised.The rate coefficients were updated to the latest recommendations of the Jet Propulsion Laboratory (JPL; S. P. Sander et al., 2011) and other recent publications (see Sect. 3.5.1).
The submodel JVAL, used to calculate photolysis rate coefficients, has been updated to the most recent version (Sander et al., 2014, see Sect. 3.5.1).

Updates of physics-related submodels
The ocean subsystem as described by Pozzer et al. (2011), consisting of the submodels MPIOM (based on the General (ocean) Circulation Model MPIOM version 1.3.0;Marsland et al., 2003), HD (Hydrological Discharge; based on the work of Hagemann and Gates, 1998;Hagemann et al., 2006), and A2O (atmosphere to ocean coupling; Pozzer et al., 2011) are now included in the new MESSy version 2.51 (see Sect. 3.11).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 (Jöckel et al., 2006), 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. 3.7.1),CLOUDOPT for the calculation of cloud optical properties, and ORBIT for Earth orbit calculations.The technical documentation with application examples will be published elsewhere (Dietmüller et al., 2016).
To increase the spectral resolution of the ultravioletvisible (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 hPa (Nissen et al., 2007;Kunze et al., 2014).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 bud-gets of the surface (surf.f90)1, 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 Roeckner et al. (2003).
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 (Roeckner et al., 2003(Roeckner et al., , 2006)), 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; Baumgaertner et al., 2013), CONVECT (convection; Tost et al., 2006b), and CLOUD.Convective tracer transport is simulated by the submodel CVTRANS (Tost, 2006).Vertical diffusion, orographic gravity waves, and large-scale advection (flux-form semi-Lagrangian; Lin and Rood, 1996) are treated in EMAC by ECHAM5.The feedback of atmospheric chemistry to the hydrological cycle is controlled by the submodel H2O (Jöckel et al., 2006).

New diagnostic submodels
New diagnostic submodels have been added: -SCALC (see also Appendix A1 for an example usage; Kern, 2013) is used for Simple CALCulations with "channel objects" (Jöckel et al., 2010) 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 Jöckel et al. (2010).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 A1) is used to calculate budgets of reactive compounds.
ing strengths are not applied homogeneously in the vertical: the boundary layer and the stratosphere-middle atmosphere above 10 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. 3.3).In addition to the RC2-base simulations (free-running hindcast and projection), which are forced by prescribed SSTs/SICs (see Sect. 3.3), we conducted a simulation with an interactively coupled ocean model (Pozzer et al., 2011).The details are described in Sect.3.11.
As listed in Table 1, 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 (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.3.7 and 3.8.Last, but not least, we performed additional sensitivity simulations (those appended by letter "a" in Table 1) for reasons outlined in Sect.3.12.
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 A. The high complexity of the applied model setups and the total number of simulation years required a considerable amount of computational resources and led to almost 2 PetaByte of primary model output (Table 2).

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; Tiedtke, 1989;Nordeng, 1994), cloud cover (submodel CLOUD; Sundqvist et al., 1989), and non-orographic gravity waves (submodel GWAVE;Hines, 1997)

84
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 m s −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 W m −2 for RC1-base-08 and RC2-base-05, respectively.
Comparison to the test simulation without interactive chemistry (0.1 W m −2 ) shows that these values are still in the range of ±0.5 W m −2 , only slightly larger than the uncertainty range from observations.Stephens et al. (2012) give an estimate for the TOA radiation balance of 0.6 (±0.4) W m −2 for the decade 2000-2010 derived from satellite observations.

Balancing the radiation budget
To achieve a balanced radiation budget at the TOA, the tuning of cloud parameters is commonly applied (e.g.Mauritsen et al., 2012).The default settings for the T42L47MA set-up of EMAC result in a TOA radiation imbalance of −3.3 W m −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 Mauritsen et al. (2012).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, possi- ble range: 0.10-0.38)and the entrainment rate for deep convection (entrpen = 0.5 × 10 −4 m −1 , default: 1.0 × 10 −4 m −1 , possible range: 0.3-3.0 × 10 −4 m −1 ).The values for the possible ranges of these two parameters have been tested by Mauritsen et al. (2012).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 W m −2 increase of the net SW radiation, and an 0.9 W m −2 increase of the OLR at TOA, which results in an almost balanced TOA radiation budget of 0.1 W m −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. 3.11 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 hPa, with full nudging weights (i.e.1.0) from 20 to 50 hPa, levelling off to 0.3 (0.2) at the upper (lower) edge of the nudging region (Giorgetta and Bengtsson, 1999).The lat-itudinal range is confined to 12.6 • S-12.6 • N, with full nudging at latitudes from 7 • S-7 • N. The relaxation timescale is set to 10 days.Although the internal generation of a QBO is a feature of the L90MA set-up of EMAC (Giorgetta et al., 2002), the zonal winds near the Equator are slightly nudged with a relaxation timescale of 58 days in all EMAC simulations at T42L90MA resolution, to get the correct phasing of the observed QBO.

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 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 (Manzini et al., 1997).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 • S and the temperatures for a polar cap average from 71.2-87.9• S (Fig. 1), 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 • 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/;Rayner et al., 2003).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. 3.1).In the global mean, the HadISST and ERA-Interim SSTs and SICs are almost identical (Fig. 11).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 ±20 % (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 (Collins et al., 2011;Martin et al., 2011).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 Jones et al. (2011).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. 11 and S10;Martin et al., 2011).The sea-ice extent has been well reproduced, remaining within 20 % of observed values during most of the year (see also Figs. 11 and S13;Martin et al., 2011), even though we find deviations of up to 80 % 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; Lean et al., 2005) have been used.Data have been prepared as input for the radiation (RAD) sub-submodel FUBRAD (Kunze et al., 2014), here applied with 55 spectral bands, and for the photolysis calculations in the submodel JVAL (see Sect. 3.5.1).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 Jones et al. (2011;see also Sect. 3.3).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).

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 R. Sander et al. (2011) and Jöckel et al. (2010), a revised version was used for the simulations.All rate coefficients were updated to the latest recommendations by JPL (S.P. Sander et al., 2011) and other recent publications.Apart from these updated rate coefficients, the product distributions were also updated for a few reactions (e.g.C 2 H 4 + O 3 ).In addition, previously neglected, chemically inert or ubiquitous products like CO 2 , H 2 O and O 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 Sander et al. (2014).Spectrally resolved time series of photon fluxes have been used consistently with the transient spectrally resolved irradiances used for the radiation calculation (see Sect. 3.4).
For the base simulations (RC1-base, RC1SD-base, and RC2-base; see Sect. 3 and Table 1), the chemical mecha-nism 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 C 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 Pöschl et al. (2000).Heterogeneous reactions of dinitrogen pentoxide (N 2 O 5 ), halogen nitrates (ClNO 3 , BrNO 3 ), and hypohalous acids (HOCl, HOBr) are also included.Since Hg chemistry is not considered in this study, all 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.3.8 and Table 1) 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 ES-CiMo_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 (Tost et al., 2006a(Tost et al., , 2007a(Tost et al., , 2010)), 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 (HNO 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 tropo-spheric values by its sub-submodel MECCA_KHET; see also Sect. 3.7.2).
For the formation of NAT, a kinetic growth NAT parameterization (Kirner et al., 2011;van den Broek et al., 2004) is used with the assumption of a necessary supercooling of 3 K (Schlager and Arnold, 1990).Hence, homogeneous NAT nucleation is also possible.The formation of STS is based on Carslaw et al. (1995), the formation of ice particles is based on the thermodynamic approach of Marti and Mauersberger (1993).For the ice nucleation 20 % H 2 O supersaturation is assumed to be necessary.A trapezoid scheme (Buchholz, 2005;Kerkweg et al., 2006a) is used for the sedimentation of PSC particles.
Heterogeneous reaction rates and their temperature dependencies on NAT are calculated according to the parameterization of Carslaw et al. (1997) based on the measurements of Hanson and Ravishankara (1993).The heterogeneous reaction rate coefficients on liquid particles are taken from Hanson and Ravishankara (1994) and Hanson et al. (1994).The uptake coefficients and reaction probabilities for ice particles are taken from S. P. Sander et al. (2011).
Here, stratospheric H 2 SO 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 • in latitude and on 70 pressure levels between 530 and 3 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 aerosolphase species.Additionally, sedimentation leads to a significant loss of aerosol particles from the atmosphere.Within the MESSy submodel DDEP (formerly named DRYDEP; see Sect.2; Kerkweg et al., 2006a), dry deposition velocities are calculated following the big leaf approach as proposed by Wesely (1989).
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 (Kerkweg et al., 2006a).
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 = 5 × 10 −07 m, σ = 2.0) and particle density (ρ = 1841.0kg m −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 (Jöckel et al., 2010).Short-lived reactive species, e.g.nitrogen oxides, nitric acid (HNO 3 ), and ozone, were initialized directly with the mixing ratios from the year 2000.During the spin-up phase (1950)(1951)(1952)(1953)(1954)(1955)(1956)(1957)(1958)(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 (CO 2 ), nitrous oxide (N 2 O), and methane (CH 4 ), but also more reactive species, such as carbon monoxide (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 (Kerkweg et al., 2006b).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 (CO 2 , N 2 O and CH 4 ), ozone depleting substances (CFCs: CFCl 3 , CF 2 Cl 2 , CH 3 CCl 3 , CCl 4 ; HCFCs: CH 3 Cl, CH 3 Br; Halons: CF 2 ClBr, CF 3 Br), H 2 , and SF 6 (Figs.E1-E6).In the RC1aero and RC1-aecl simulations COS has been prescribed as well (Brühl et al., 2012, see Fig. E7).
For all species except 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 SF 6 were taken from Meinshausen et al. ( 2011) and extended with the RCP6.0 scenario as proposed by Eyring et al. (2013b).Mixing ratios of ozone depleting substances were 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 H 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 RC1SDbase-10a (see Table 1 and Sect.3.12.2) we used the RC2 setup 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. 3.8) 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 us-ing tropospheric values below 500 hPa, stratospheric values above 300 hPa, and linearly interpolated values between 500 and 300 hPa.
For the troposphere the Tanre climatology (Tanre et al., 1984), 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 lowresolution 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 H 2 SO 4 mixing ratio (monthly zonal averages, 5 • latitude resolution; see Sect.3.5.3).The corresponding values were inter-/extrapolated to the radiation bands of the EMAC model and provided with a 500 m grid spacing in the vertical.As for H 2 SO 4 , data from the year 1960 have been used for the spinup 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 N 2 O 5 into HNO 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 Righi et al. (2013).This simulation was performed with EMAC (MESSy version 1.4) coupled to the aerosol submodel MADE.It covers a period of 10 years (1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(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 (Lamarque et al., 2010).
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. 3.5.3).

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 RC1aecl simulations.To determine the physical and chemical properties of the atmospheric aerosol, the EMAC aerosol submodel GMXE (Pringle et al., 2010;Tost and Pringle, 2012) 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 (Fountoukis and Nenes, 2007), 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 Jöckel (2006) in the new submodel IM-PORT (Kerkweg and Jöckel, 2015), vertically redistributed via OFFEMIS (Kerkweg et al., 2006b), 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 (Tost et al., 2006a;Kerkweg et al., 2006a).The RC1-aero and RC1-aecl simulations consider the effect of phase partitioning of H 2 SO 4 , HNO 3 , HCl, and NH 3 , as well as the interaction of these compounds with primary aerosol species such as 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 (Pozzer et al., 2012;Klingmüller et al., 2014), replacing the climatology in the troposphere.Note that the extensions of Klingmüller et al. (2014) 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 Chang et al. (2014).The activated aerosol number is subsequently used to determine prognostic cloud droplet number concentrations in the two moment cloud microphysics scheme of Lohmann and Ferrachat (2010).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 2).

Prescribed emissions
Anthropogenic emissions are incorporated as prescribed emission fluxes following the CCMI recommendations (Eyring et al., 2013b).Two data sets are considered: one data set is the MACCity3 (Granier et al., 2011;Diehl et al., 2012;Lamarque et al., 2010) 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 (Lamarque et al., 2010) and RCP 6.0 data (Fujino et al., 2006).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 Pozzer et al. (2009).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 nonanthropogenic emissions such as NMHCs of biogenic origin, terrestrial dimethyl sulfide (DMS), volcanic SO 2 , NH 3 , halocarbons, CH 3 I, and OC from secondary organic aerosol (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.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 1 and Sect.3.12.2) we used the RC2 set-up for the years 2012 and later.

Biogenic emissions
The emissions of NO x from soil and isoprene (C 5 H 8 ) from biogenic sources are calculated on-line using the submodel ONEMIS.Besides prescribed fields, like the distribution of cultivation and agriculture, the calculated NO x emissions depend mainly on the soil temperature and the soil wetness.The algorithm used is based on Yienger and Levy (1995) as described by Ganzeveld et al. (2002).Estimates for the soil biogenic emissions of NO are highly uncertain.Vinken et al. (2014) compared results of different models and satellite observation-based estimates.Most of the annual totals are in the range of 4 to 15 Tg (N) a −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. 2. In all simulations the emission totals are between 5.4 and 6.3 Tg a −1 until the year 2010 and thus at the lower end of the estimated range from Vinken et al. (2014).While the simulations of the RC2 series are at the lower limit of this range, the simulations RC1SDbase-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 Tg (N) a −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 Guenther et al. (1995) and implemented according to Ganzeveld et al. (2002).For the biogenic isoprene emissions (see Fig. 3) we simulate a range between 430 and 550 Tg (C) a −1 .Similar as for the biogenic 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 Tg (C) a −1 until 2100.However, Guenther et al. (1995) do not consider the CO 2 inhibition effect on isoprene emissions, implying that the future isoprene emissions (as projected in the RC2 simulations) may be too strong.Guenther et al. (2006) estimated the annual total emissions of isoprene from biogenic origin to be 440 to 660 Tg (C) a −1 (see also Arneth et al., 2008, and 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 Jöckel et al., 2006).A more detailed discussion on the isoprene scaling factor is provided by Pozzer et al. (2007).

Lightning NO x
The NO x emissions from lightning activity are calculated online using the MESSy submodel LNOX (Tost et al., 2007b).Here, we apply the parameterization by Grewe et al. (2001), 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 NO x emissions by lightning differ between the simulations.
Estimates for the annual total emissions of NO x from lightning are in the range of 2-8 Tg (N) a −1 (see Schumann and Huntrieser, 2007).In Fig. 4 we see that most of the simulations are within this range.The results cluster around values between 4 and 5 Tg (N) a −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. 4.1).Despite the large difference of the annual total lightning 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).

Ocean-to-atmosphere fluxes
Ocean-to-atmosphere fluxes of DMS, C 5 H 8 , and methanol (CH 3 OH) are calculated by the AIRSEA submodel (Pozzer et al., 2006).The ocean-to-atmosphere flux of a chemical species is calculated from its concentrations in the upper- most ocean layer and the lowermost atmosphere layer following the two-layer model by Liss and Slater (1974).We use the parameterization of Wanninkhof (1992) for the water side exchange velocity in this study.The parameters needed for the simulated species in AIRSEA correspond to the values suggested by Pozzer et al. (2006).Ocean salinity is taken from the monthly climatology of the World Ocean Atlas 2001 (Boyer et al., 2002).Oceanic concentrations of DMS are a monthly climatology from Lana et al. (2011).Isoprene concentrations in the ocean are calculated using the parameterization of Broadgate et al. (1997), relating isoprene and chlorophyll concentrations, here with chlorophyll prescribed as a monthly climatology from the World Ocean Atlas 2001 (Conkright et al., 2002).For methanol, the atmosphere-toocean flux is calculated assuming a constant undersaturation of surface water methanol with respect to lowermost atmospheric concentrations (Singh et al., 2003).The constant saturation coefficient for methanol in surface water is 0.94.
The annual global ocean-to-atmosphere fluxes of DMS are shown in Fig. 5; total oceanic emissions range between 25.5 and 31.5 Tg (S) a −1 until 2010, and up to 32.5 Tg (S) a −1 in the projection simulations.These values are within the range reported in the literature of 15 to 54 Tg (S) a −1 (Kettle and Andreae, 2000, and references therein).The trend in DMS emissions from the ocean follows the trend of SSTs (Fig. 11).Therefore, compared to the -base-simulations, the increase in emissions is lower in the coupled ocean simulation RC2oce-01, as the latter shows a lower increase of SSTs.
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 (Erickson and Hernandez, 2013) to  (Palmer and Shaw, 2005), much lower than previous estimates (Bonsang et al., 1992).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 Tg (C) a −1 (Fig. 6), 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 (Conkright et al., 2002) used to estimate the isoprene concentration in the surface water.The trend of isoprene emissions follows the trend of SSTs (Fig. 11).
For methanol, the net oceanic sink is uncertain, ranging from 0.1 to 21 Tg (C) a −1 (Heikes et al., 2002;Galbally and Kirstine, 2002;Jacob et al., 2005), with a recent study by Millet et al. (2008) resulting in an net oceanic methanol sink of 6 Tg (C) a −1 .The simulated oceanic uptake of methanol ranges from 1.2 to 1.85 Tg (C) a −1 (Fig. 7).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 % (Jacob et al., 2005), whereas we assumed 94 %.Furthermore, the net oceanic sink for methanol is still uncertain and only one publication reporting measurements of methanol concentrations in seawater exists (Williams et al., 2004).
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 Wanninkhof (1992) 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 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 Yang et al. (2005), we assume that 50 % of the bromide is released to the gas phase (leading to an additional factor of 0.5).The resulting time series of oceanic Br emissions is shown in Fig. 8.The absence of a trend in sea-salt emissions indicates that the 10 m wind speed over the ocean does not change significantly on global and annual average.This is consistent with the study by Pryor et al. (2006), revealing that changes in wind speeds from general circulation model (GCM) projections are small (≤ 15 %).
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 Burrows et al. (2013).The result is shown in Fig. 9.The oceanic organic carbon source is highly uncertain, with recent estimates of 8 Tg (C) a −1 (Spracklen et al., 2008) and 75 Tg (C) a −1 (Roelofs, 2008).Our simulated annual total oceanic emissions of POC range between 17.9 and 19.5 Tg (C) a −1 , which is close to the first estimate of 14 Tg (C) a −1 for oceanic POC emissions by Duce (1978).

Dust emissions
Dust emissions are calculated on-line in the RC1-aero and RC1-aecl simulations.We use the dust emission scheme of Tegen et al. (see Figs. E58 and E59 for input parameters;2002), which was implemented in the submodel ONEMIS by Gläser et al. (2012).The wind stress threshold for dust emissions is corrected by a factor of 0.86 in accordance with Tegen et al. (2004).The on-line calculated total global dust emissions are in the range between 600 and 1400 Tg a −1 (Fig. 10), of which 1.5 to 2 % 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 Tg a −1 (Tegen et al., 2004), dust emissions of the RC1-aecl simulations are within the suggested range.
As shown in Fig. 10, 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. 14) 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.
The simulations with prescribed aerosol use the Tanre climatology (Sect.3.7.1),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 atmosphereocean 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 • × 3 • , 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 Kern (2013).In addition to the parameters described in Sect.3.2, 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 (CO 2 , CH 4 , N 2 O, 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 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 preindustrial conditions, GHGs were prescribed representative for the year 1750, i.e. without any trend.The applied values are 0.28×10 −3 , 0.72×10 −6 , and 0.27×10 −6 mol mol −1 for CO 2 , CH 4 , and N 2 O, respectively.The mixing ratios of the CFCs were set to zero.The remaining radiation imbalance in SP-oce-01 averaged over the last 30 years is 1.375 W m −2 with a globally averaged surface temperature of 288.4 K.This imbalance seems to be quite high; however, the corresponding optimized parameter set is the result of a multitude of sensitivity studies (Kern, 2013), yielding the optimum results in terms of climate and hydrological cycle for both, preindustrial and industrial conditions.At the end of the spin-up period of SP-oce-01, no statistically significant trends (90 % confidence level) of the surface temperature over the last 60 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 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, O 3 has been prescribed from a monthly climatology, which has been derived from the years 1962-1972 of a previous test simulation.Austin et al. (2012) 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 Nowack et al. (2015).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 RC2base-08.
Figure 11 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 (IPCC, 2013, Fig. AI.4), which show an ensemble median (average) increase of 1.7 (2.1) K for the same scenario (RCP6.0).Nevertheless, this increase (about 1.5 K, Fig. 11, upper panel) is below what has been projected by the HadGEM-ES model (about 2.5 K), which also implies a slower decrease of the global sea-ice coverage (Fig. 11, lower panel).Despite this smaller increase in temperature, the results are still in line with the CMIP5 simulations (Collins et al., 2013;IPCC, 2013), where the projected increase is between 1.3 and 2.7 K.Moreover, the HadGEM model is at the high end of the CMIP5 models in terms of climate sensitivity (Andrews et al., 2012).

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 (Pipitone and Easterbrook, 2012).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 2) 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.3.12.3),we decided for consistent time series, but performed additional sensitivity simulations (indicated by suffix a; see Table 1) 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.3.7.1)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 multiannual monthly average radiative heating rates  from simulations RC1-base-08, RC1-base-08a, RC1-base-07, and RC1-base-07a.Above 100 hPa the range of differences between corresponding pairs is smaller than the interannual 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.3.5.3),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.3.12.1)corrected.The corrections cause an increase of the total ozone column of up to 2 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 RC1aero-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.Grewe et al., 2012) and without a detailed analysis of statistical implications on the evaluation of the simulation data (e.g.Grewe and Sausen, 2009).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.

Temperature profiles
For the assessment of the simulated air temperatures, we compare the simulation results with ERA-Interim data (Dee et al., 2011).We choose the period 2000-2010 for the evaluation, as it is covered by most of the simulations.Multiannual climatologies of zonally averaged temperature profiles for the period 2000-2010 are shown in Fig. 12. Figures 12 and S18 were produced with the ESMValTool (Righi et al., 2015;Eyring et al., 2015).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 freerunning 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 (Brinkop et al., 2015).
The simulations in the first category perform best, as expected, compared to ERA-Interim with minor differences of mostly less than ±1 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 RC1SDbase-09 simulation is at around 70 hPa with more than ±4 K. Between 100 and 150 hPa this simulation shows a small bias of ±1 to ±1.5 K.In the RC1SD-base-10 simulation, with a higher vertical resolution, the maximum of the cold bias is located lower, at about 200 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 RC1base-07 in Fig. 12).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 K) around the tropopause and show a warm bias (of up to 4 K) in the Southern Hemisphere (SH) above 100 hPa, as also analysed from previous EMAC simulations by Righi et al. (2015).
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 • S-30 • N) tropopause at around 250 hPa.The freerunning simulations show a more pronounced seasonal cycle of the bias than the nudged simulations.The warm bias above the tropopause around 60-30 • S, which is also seen in the annual climatology, is strongest in SH winter (more than 4 K).Around this time there is also a strong cold bias (of up to 4 K) in the tropics above the tropopause.

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/;Dee et al., 2011) 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;Adler et al., 2003) data are available as monthly means for each month since January 1979 at 2.5 • × 2.5 • horizontal resolution.This precipitation data set combines satellite data with rain gauge measurements, resulting in a data set with global coverage.As Dee et al. (2011), 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. 13.
All simulations, as well as ERA-Interim, overestimate precipitation in the tropics by 0.14 up to 0.84 mm day −1 , depending on the simulation (see Fig. 13).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 mm day −1 ) in the NH and minima (about 2 mm day −1 ) of precipitation at approximately 25 • N and 25 • 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 mm day −1 ) are reproduced by ERA-Interim and all EMAC simulations (about 3 to 3.75 and 2.25 to 2.75 mm day −1 ).Apart from the tropics, the largest deviations can be found in the SH subtropics.Horizontal distributions of average precipitation deviations (in mm day −1 ) throughout the years 1990-2009 are shown in Fig. 14.
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. 13) 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. 13) 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 ±3 mm day −1 ) are found in the tropics, where rainfall is the strongest.
Overall the magnitude of deviations is comparable to uncertainties that arise when changing convection parameterizations within ECHAM5 (Tost et al., 2006b) and have also been found in the analysis of the ECHAM5 model (Hagemann et al., 2006).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 (Dai, 2006).

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. 15 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 Tg a −1 , the deposition fluxes increase up to a maximum of 850 to 880 Tg a −1 around the year 2000.This is well within the range of other chemistry climate model results.For instance, Young et al. (2013) reported ozone deposition fluxes in the range between 687 and 1360 Tg a −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 nmol mol −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 NO x and corresponding ozone production (see Fig. 4 and Sect.4.7).The wet deposition fluxes of ozone exhibit the same features, but are 6 to 7 orders of magnitude smaller (see Fig. S21).
For the NO x wet deposition fluxes the difference between L47MA and L90MA amounts to 20 % (see Fig. S22).In contrast to ozone, the NO x dry deposition fluxes peak earlier (around 1990) and decrease slowly afterwards (see Fig. S22).
The deposition fluxes for nitrate (HNO 3 + N 2 O 5 + NO − 3 ) and sulfate (H 2 SO 4 + HSO − 4 + SO 2− 4 ) are shown in Figs.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.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 sim-ulations 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 coarsemode aerosol radii: 0.5 µm (prescribed); simulated: 0.759 µm (RC1-aero-06), 0.751 µm (RC1-aero-07), 0.793 µm (RC1aecl-01), 0.782 µ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 processesthe 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 (Kerkweg et al., 2006a).
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. 17 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. 17, 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 % smaller, the wet deposition fluxes up to 15 % smaller inaero-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 % 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 (Lamarque et al., 2013) for the year 2000 sulfur emissions (especially DMS) are slightly enhanced, whereas NO x and NH 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 Lamarque et al. (2013), with, e.g., 0.64 for oxidized sulfur compared to 0.72.Similar differences are found for NH x and 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 (Tost et al., 2007a).

Tropospheric oxidation capacity
The hydroxyl radical (OH) is the most important cleansing oxidant in the atmosphere.However, variations in OH concentration are still difficult to quantify (Manning et al., 2005;Patra et al., 2014;Ghosh et al., 2015).Among others, it reacts with CH 4 largely determining its atmospheric lifetime.
Here, we present the simulated OH lifetime of atmospheric CH 4 at time t as a measure for the oxidising power of the atmosphere, calculated according to Jöckel et al. (2006) as 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 ± 0.6 years in the period of 2000-2004 (mean and standard deviation for all simulations covering this period).In comparison to other methane lifetime (versus OH) estimates, as by Voulgarakis et al. (2013), Jöckel et al. (2006), von Kuhlmann (2001), Hein et al. (1997), andRighi et al. (2015) (9.8 ± 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.Prather et al., 2012;Prinn et al., 2005, with 11.2 ± 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 (Lawrence et al., 2001), whereas varying included vertical layers due to different tropopause heights have a minor impact (see also O'Connor et al., 2014).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 CH 4 by OH and also has a positive impact on the 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 OH mixing ratio in RC1SD-base-07, mainly in the tropical upper tropopause, caused by a likewise smaller NO x production from lightning activity (see Fig. 4).The longest lifetime is simulated in RC1-aecl-02 (burlywood), which shows the largest cold biases of all simulations.
Moreover, the variations in production of OH by tropospheric ozone depend on various factors (Hofzumahaus et al., 1992).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 Brenninkmeijer et al. (2007), 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 ((CH 3 ) 2 CO).CARIBIC flies to destinations worldwide; thus in order to obtain meaningful climatologies with good data coverage, latitude was limited to 35 • N < latitude < 60 • N. Data were also limited by including only values where pressure p < 280 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 • grid resolution) as distance to the 3.5 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.
Results of the data prepared in this way are presented in Fig. 19 for O 3 , Fig. 20 for CO, Fig. 21 for CH 4 , and Fig. 22 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: RC1SDbase-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 O 3 (Fig. 19), the model results compare well with measurements (see Zahn et al., 2012, for instrument specifics).The seasonal cycle is reproduced, but the lower values in the troposphere are generally overestimated by up to 40 % by the model.In the stratosphere, differences are smaller, as the model underestimates measurements by 5 %, reaching 30 % 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 O 3 (shown in Supplement).When nudging global mean temperature, upper tropospheric O 3 levels increase by about 5 %, while stratospheric levels stay constant.
Increasing the vertical resolution on the other hand seems to decrease the differences between the seasons, reducing O 3 levels in spring and summer by about 10 % in the simulations with 90 levels and increasing them by up to 5 % in the Contrary to O 3 , it is the increasing vertical resolution that produces differences between stratosphere and troposphere for CO, tropospheric values increasing by 5 % for the simulations with 90 levels, and being reduced by up to 10 % in the stratosphere.
Results for CH 4 (Fig. 21; see Dyroff et al., 2014, for instrument specifics) are similar.Relative differences in general diff.
[%] 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. 4.4).CH 4 is noticeably affected by the road traffic emissions, decreasing the mean values by about 1 % randomly, not showing any seasonality or stratification (shown in Supplement).This is consistent with a decreased CH 4 OH lifetime (Fig. 18).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 % in the troposphere and down to −0.5 % in the stratosphere for model simulations with an increased Acetone meas, @35<lat<60 @wholeday @180<p<280 vertical resolution (see Sect. 4.4).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 %.Acetone (Fig. 22) shows the largest deviations, when comparing model results with measurements (for instrument specifics, see Neumaier et al., 2014).Tropospheric values are underestimated by about 60 to 100 % 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 % in the stratosphere during late winter and spring (January to May).Increasing the vertical resolution also increases mean acetone values by 5 % in all seasons and at all heights, except for summer and autumn (May to November) stratospheric values, which are reduced by up to 10 %.These results are in line with the previous evaluation of Pozzer et al. (2007, Table 8), who observed an underestimation of ∼ 50-35 % for the same biogenic emissions as used here (55 Tg a −1 ), which account for ∼ 85 % of the total emissions.These results indicate the need for considering oceanic emissions and photochemical production of CH 3 COCH 3 from monoterpenes, methylbutenol, and higher iso-alkanes (Jacob et al., 2002;Khan et al., 2015), which were not included.For example Pozzer et al. (2010) showed that i-butane and i-pentane are responsible for the photochemical production of 4.3 and 5.8 Tg a −1 of CH 3 COCH 3 , respectively.

Measurement
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 %.Different strengths of this influence can be noted for most species between the stratosphere and the troposphere, for O 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 O 3 , leading to relative differences of up to 10 % for selected months.Overall, the model produces realistic values for O 3 , but underestimates CH 4 , 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, quasihorizontal mixing processes contribute to the stratospheric transport (Plumb, 2002).The RC is often expressed in terms of the tropical upward mass flux (F 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 (Holton, 1990).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; Hall and Plumb, 1994).It is obtained from linearly increasing surface emissions of an inert tracer, e.g.sulfur hexafluoride (SF 6 ), which are mostly used from observations.The AoA contains both components of the BDC, RC and mixing processes.
Results on F trop and mean AoA are shown in the Sect.4.6.1 and 4.6.2.For evaluation of the model results on F trop , the RC is calculated from ERA-Interim reanalysis.Moreover, data of atmospheric SF 6 concentrations from July 2002 to April 2012 are available from the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS; Fischer et al., 2008) instrument for the calculation of mean AoA.Here, we use mean AoA from the MIPAS-ENVISAT (Environmental Satellite) level 1b spectra version 5 (Haenel et al., 2015), which is an update of Stiller et al. (2012).

The tropical upward mass flux
Figure 23 shows the time evolution of the annual mean tropical upward mass flux from 1960 to 2100 from the simulations in the lower (around 70 hPa, top), middle (around 10 hPa, middle), and upper (around 1 hPa, bottom) stratosphere.ERA-Interim data are included from 1979 to 2013.All simulations show an increase in F trop with rising atmospheric mixing ratios of GHGs, which is well known from earlier CCM simulations (e.g.Butchart et al., 2010;Oberländer et al., 2013).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 hPa), which is not captured by the model simulations.As shown by Seviour et al. (2012), this negative trend contains large uncertainties, and does not occur in other reanalysis systems or when using different estimates of upwelling in ERA-Interim (Abalos et al., 2015).The nudged simulations, especially the simulations RC1SDbase-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.Garny et al., 2011;Oberländer et al., 2013).The simulation RC2-oce-01 with interactively coupled ocean (Sect. 3.11) 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 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 hPa (Sect.3.1), the RC1SD-simulations do not perform better than the other simulations.

Mean age of stratospheric air
Figure 24 shows the mean age of stratospheric air from the simulations (colours) and MIPAS data (black) averaged from 2002 to 2011 at around 50 hPa (20 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 MI-PAS observations, the L90 results fit well with a mean AoA of around 3 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 SF 6 measurements by MIPAS is affected by the mesospheric sink of SF 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 SF 6 of up to 2 years (Haenel et al., 2015).A lower mean AoA and therefore a faster mass transport from model simulations compared to observations is well known from other CCMs (SPARC, 2010).The simulation with coupled ocean lies in between the vertically lower and higher-resolved simulations and thereby fits better with observations from MI-PAS 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. 25.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 trop (Sect.4.6.1),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.
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. 25, 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 lowerresolved simulations.In contrast to the findings for the RC (Sect.4.6.1), the vertical resolution is much more important for the mean AoA than different SST data sets.
In the middle and upper stratosphere (Fig. 25, 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 SF 6 tracer for the observations.SF 6 is inert in most parts of the stratosphere, but photochemically destroyed in the mesosphere.The decent of SF 6 -poor air from the mesosphere into the stratosphere raises the calculated mean AoA (Stiller et al., 2012).This might explain some of the deviations between EMAC and MIPAS data in the upper stratosphere.
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; Bodeker et al., 2005) for total column ozone (TCO), AURA Microwave Limb Sounder/Ozone Monitoring Instrument (MLS/OMI; Ziemke et al., 2011) for tropospheric and stratospheric partial column ozone (TPCO, SPCO), and the ozone sonde data set described by Tilmes et al. (2012) 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 RC1aero-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 Righi et al. (2015), and also share some routines with the ESMValTool (Eyring et al., 2015).
Figure 26 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 cal- culations 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 % 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).
Figure 27 shows the climatological annual cycle of TCO for the BSTCO observations (Fig. 27a) and the highresolution (L90MA) RC1 (Fig. 27b), RC2 (Fig. 27c) and RC1SD (Fig. 27d) simulations as well as the simulation with coupled ocean (RC2-oce-01, Fig. 27e).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. 27f-i) show that in SH spring the positive bias related to the underestimation of the ozone hole is the largest in the free-running RC1base-07 (Fig. 27f) and RC2-oce-01 (Fig. 27i) simulations, whereas the smallest bias is found in the nudged simulations (Fig. 27h).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. 12.We found that (planetary scale) wave fluxes in SH mid-latitudes are overestimated, and thus the polar vortex is too disturbed.The the 1990s (e.g.World Meteorological Organisation, 2014), 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 extratropical TCO with a return to 1980 levels until 2100.This positive trend is consistent with earlier CCM studies (e.g.SPARC, 2010;World Meteorological Organisation, 2014).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. 4.6.1)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. 29).Regarding the freerunning 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. 30.It is positive in all regions and simulations (see also Righi et al., 2015).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 K (see Sect. 4.1), the temperature-dependent reaction kinetics is altered, pushing the chemistry into a new equilibrium state (e.g.Rasmussen et al., 2013).
(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 hPa) of RC1SD-base-07 (with T nudging) compared to RC1SD-base-10 (without T nudging) reduces the TPCO polewards of 40 • latitude by up to −0.4 DU (Fig. S38).Around 30 • N/S, the corresponding shift is up to −5 hPa increasing the TPCO by up to 1.5 DU.The resulting global effect of this geometrical tropopause shift (analysed for the 11-year period 2000 to 2010) is −0.1 DU.This effect is smaller, however, than the influence of the tropopause definition: using the online diagnosed tropopause (WMO definition equatorward of 30 • latitude; iso-surface of 3.5 PVU potential vorticity poleward of 30  This is based on their finding that the average ozone production efficiency by NO x from lightning is 5 times higher compared to other NO x sources (except for aviation), and that 70 Tg (N) a −1 on average account for about 26.5 DU ozone (i.e.0.4 DU (Tg (N) a −1 ) −1 on average for all sources).This dominating effect of the modified lightning 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 NO x production peaks (Fig. S40, Appendix A5).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 (Meinshausen et al., 2011).
Vertically resolved ozone is compared to ozone sonde data (Tilmes et al., 2012).This analysis has been adapted from Righi et al. (2015), also sharing some routines with the ES-MValTool (Eyring et al., 2015).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 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 % level), but each t test is based on 12 data pairs only.Figure 31 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 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 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 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  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 km altitude (in the MA), respectively.One simulation  in T42L47MA was performed with an interactively coupled ocean model and set up based on an extensive spinup 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 PetaByte comprising data set.
First analyses presented here focus mainly on an intercomparison between the different simulations from a global perspective.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 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 • S) above 140 hPa and in the SH polar region between 140 and 30 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 DU increasing from about 5 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 trop-ics 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 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 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 DU in tropospheric total ozone column per Tg (N) a −1 lightning NO x production.The reduced lightning NO x production with global mean temperature nudging is partly compensated by an increased NO x release from soil, due to increased soil temperatures.Whereas 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 (DMS, C 5 H 8 , Br), in contrast, are reduced by nudging, and even more reduced, if global mean temperature is nudged.The oceanic uptake of CH 3 OH is likewise reduced with (global mean temperature) nudging.Deposition fluxes through scavenging and dry deposition of ozone, nitrate, 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 NO x this is consistent to the likewise-reduced lightning NO x production and the reduced O 3 burden.
The increased tropospheric temperature through the global mean temperature nudging slightly reduces the tropospheric methane lifetime towards OH, this effect is larger in simulations with the coarser vertical resolution (L47MA).Overall, i.e. averaged over all simulations, the simulated OH lifetime of methane is 8.0 ± 0.6 years in the period 2000-2004.This is potentially too short, indicated by the consistently underestimated CH 4 mixing ratios in the NH UTLS, i.e. compared to observations at CARIBIC flight levels.Note that CH 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.RC1SDbase-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 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.This submodel analyses the contribution of different production cycles (e.g.associated with different source gases) to selected tracers (type-2 family, see Jöckel et al., 2008).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 d dt 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 Br y and Cl y are defined in the CTRL_FAMILY namelist of submodel TRACER, submodel SCALC (Kern, 2013) is used to sum the loss rates LossBr and LossCl from their individual tendencies calculated by SCAV, and in the CPL namelist of TBUD-GET the corresponding diagnostic tracers BrS, BrL, ClS, and ClL for short (S)-and long (L)-lived halogen compounds, respectively, are defined for integration according to Eq. (A1).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): ProdLBr, ProdSBr, ProdLCl, ProdLBr.Note that these diagnostic tracers contain the accumulated (over time) rates.

A3 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 (Schumann, 1996).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 Burkhardt et al. (2008) and were output (channel "contrail_gp") as 10-hourly global snapshots.

A4 Specific sampling of model data
With the submodel SCOUT (Selectable Column OUTput; Jöckel et al., 2010, 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, Jöckel et al., 2010, 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; Jöckel et al., 2010, Sect. 5.1) isentropes of 340, 380, and 420 K and iso-surfaces of 2, 3, and 4 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 • latitude, and as 3.5 PVU iso-surface poleward.All corresponding surfaces and maps are output as 10-hourly snapshots.
With the submodel S4D (sampling in 4 dimensions; Jöckel et al., 2010, 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).MECCA Stratospheric ozone tracer; Newtonian relaxation (TNUDGE) j towards O 3 in the stratosphere; destroyed in the troposphere as O 3 ; the corresponding loss rate tracer LO3(s) is a qualitative measure for the troposphere to stratosphere exchange of ozone (Roelofs and Lelieveld, 1997;Jöckel et al., 2006) a The Supplement related to this article is available online at doi:10.5194/gmd-9-1153-2016-supplement.

Figure 2 .
Figure 2. Annual total emissions of on-line calculated biogenic/soil NO x emissions in Tg (N) a −1 .

Figure 3 .
Figure 3. Annual total emissions of on-line calculated biogenic isoprene emissions in Tg (C) a −1 .

Figure 4 .
Figure 4. Annual total emissions of lightning NO x in Tg (N) a −1 .

Figure 5 .
Figure 5. Simulated annual total emissions of DMS (in Tg (S) a −1 ) from the ocean.

Figure 6 .
Figure 6.Simulated annual total emissions of isoprene (in Tg (C) a −1 ) from the ocean.

Figure 7 .Figure 8 .
Figure 7. Simulated annual total flux of CH 3 OH (in Tg (C) a −1 ) between atmosphere and ocean.The negative sign indicates a net uptake by the ocean.

Figure 9 .
Figure 9. Simulated annual total emissions of POC (in Tg (C) a −1 ) from the ocean.

Figure 12 .
Figure 12.Upper left panel: climatology of annual average total dry air temperature of ERA-Interim in K.The data were monthly and zonally averaged for the period 2000-2010.Other panels: dry air temperature differences (in K) of the simulations compared to ERA-Interim data.The differences, unless grey shaded, are significant on a 95 % confidence level according to a two-sided Welch's test. 16 Figure 14.(a-f) Mean precipitation differences (mm day −1 ) for the 20-year period 1990-2009.The differences show the simulation results minus GPCP data.

Figure 15 .
Figure 15.Globally integrated, annual dry deposition flux of ozone in Tg a −1 .
where M b CH 4 being the mass of methane, k b CH 4 +OH the reaction rate of the reaction CH 4 + OH, c b air (t) the concentration of air, and OH b (t) the mole fraction of OH in the grid box b ∈ 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 OH and on the other hand by temperature, since the reaction rate of CH 4 + OH depends on temperature.Figure 18 shows the simulated tropospheric methane lifetimes (towards 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 RC2base-05 simulations (orchid, solid, and dashed lines, respectively) predict a similar decreasing behaviour for the future.

Figure 18 .
Figure 18.Annually averaged methane lifetime with respect to OH of the different simulations calculated with Eq. (1).

Figure 19 .
Figure 19.Comparison of O 3 climatologies (35-60 • 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 km) relative to the tropopause.

Figure 22 .
Figure 22.As Fig. 19, but for acetone.The maximum of measured acetone values in summer is 1.5 nmol mol −1 .

Figure 26 .
Figure 26.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 28 .
Figure 28.Time series of annual mean total column ozone (in DU) averaged over different latitude bands.BSTCO observations are shown in black.Note the different scales for the different regions.

Figure 29 .
Figure29.Left column: climatological annual mean of observed total column ozone (BSTCO, top) as well as (middle) and tropospheric (bottom) partial column ozone (AURA MLS/OMI), all in 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 % confidence level are coloured.

Figure 30 .
Figure 30.Same as Fig. 26, 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.

Table 1 .
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.

Table 1 .
Overview of ESCiMo Reference simulations.The CCMI notation Ref-C1/2 is abbreviated with RC1/2, respectively, 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 setup).An appended letter "a" indicates a sensitivity study (see remarks and Section 3.12).For detailed explanations, see text.The line colours refer to the figures below.

Table 2 .
Used computational resources and size of simulated data a Including restart files (every 3 months and QTIMER triggered at the end of the scheduler wall-clock limit); b on 4 nodes of an IBM Power6 in SMTP mode (i.e. with 64 tasks/node); c temporal overlap of 06 with 07: January 1990-August 1998; d temporal overlap of 01 with 02: January 1965-October 1972.
where it is set to 0.96; for the zonal mean zonal winds at 60 • S in m s −1 (top); for the zonal mean temperatures averaged from 71.2 to 87.9 • S in K (bottom).The shading indicates the differences that are significant at the 95 % level, estimated with a Student's t test.
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Figure 1.Differences of the daily mean climatology (19 years) between a T42L47MA EMAC simulation with the gravity wave parameter rmscon set to 0.92 and a simulation

www.geosci-model-dev.net/9/1153/2016/ Geosci. Model Dev., 9, 1153-1200, 2016 1188 P. Jöckel et al.: ESCiMo with MESSy v. 2.51
(Tilmes et al., 2012)terns of all simulations are Same as Fig.26, but for the ozone volume mixing ratio at 6 selected pressure levels (700, 500, 250, 100, 50, 10 hPa).The ozone sonde data set(Tilmes et al., 2012)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.

Table A1 .
Additional diagnostic tracers.Acronyms in parentheses denote the submodels simulating the indicated processes.Corresponding data files are listed in Table E1.Newtonian relaxation (TNUDGE) c towards time series based on observations (see Fig. E6) AOA PTRAC Age of air tracer; Newtonian relaxation (TNUDGE) c at lowest model layer towards linearly in time increasing mixing ratio SF6_AOA PTRAC Age of air tracer; Newtonian relaxation (TNUDGE) c at lowest model layer towards a latitude dependent, linearly in time increasing mixing ratio SF6_AOAc PTRAC Age of air tracer; Newtonian relaxation (TNUDGE) c at lowest model layer towards a linearly in time increasing mixing ratio SF6_CCMI e PTRAC Emissions (OFFEMIS) according to EDGAR v4.2 database (see Fig. E55) SO 2 t f,g PTRAC anthropogenic emissions (OFFEMIS) as SO 2 , wet removal (SCAV) as SO 2 NH_05 d (Eyring et al., 2013b; Jöckel et al. (2010; Tracer Release EXPeriment; Sect.6.3);crelaxationtime constant τ = 3 h; d according to CCMI(Eyring et al., 2013b; Sect.4.2); e dissenting from d without interpolation of emission time series to monthly values; f dissenting from d with transient anthropogenic emissions; g see ESCiMo_SCAV_mechanism.pdf in the Supplement; h dissenting from d with seasonal cycle of emission; i see ESCiMo_MECCA_mechanism.pdf in the Supplement; j relaxation time constant is equal to model time step length.