the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
ATTILA 4.0: Lagrangian advective and convective transport of passive tracers within the ECHAM5/MESSy (2.53.0) chemistry–climate model
Patrick Jöckel
We have extended ATTILA (Atmospheric Tracer Transport in a LAgrangian model), a Lagrangian tracer transport scheme, which is online coupled to the global ECHAM/MESSy Atmospheric Chemistry (EMAC) model, with a combination of newly developed and modified physical routines and new diagnostic and infrastructure submodels. The new physical routines comprise a parameterisation for Lagrangian convection, a formulation of diabatic vertical velocity, and the new gridpoint submodel LGTMIX to calculate the mixing of compounds in Lagrangian representation. The new infrastructure routines simplify the transformation between gridpoint (GP) and Lagrangian (LG) space in a parallel computing environment. The new submodel LGVFLUX is a useful diagnostic tool to calculate online vertical mass fluxes through horizontal surfaces. The submodel DRADON was extended to account for emissions and changes of ^{222}Rn on Lagrangian parcels. To evaluate the new physical routines, two simulations in freerunning mode with prescribed sea surface temperatures were performed with EMAC–ATTILA in T42L47MA resolution from 1950 to 2010. The results show an improvement of the tracer transport into and within the stratosphere when the diabatic vertical velocity is used for vertical advection in ATTILA instead of the standard kinematic vertical velocity. In particular, the ageofair distribution is more in accordance with observations. The global tropospheric distribution of ^{222}Rn, however, is simulated in agreement with available observations and with the results from EMAC in grid space for both Lagrangian systems. Additional sensitivity studies reveal an effect of interparcel mixing on the age of air in the tropopause region and the stratosphere, but there is no significant effect for the troposphere.
 Article
(5179 KB) 
Supplement
(291 KB)  BibTeX
 EndNote
Due to the increasing demand for including interactive tracers in climate simulations it is becoming necessary to use global models which meet the needs of a fast and exact tracer transport scheme. Commonly used methods to describe largescale transport in a general circulation model of the atmosphere follow the Eulerian method. The Lagrangian (LG) method (i.e. from the perspective of a fluid particle or parcel) is more frequently used offline for trajectory studies in particle models like the global 3D chemistry transport model (Collins et al., 2002), FLEXPART (Stohl et al., 1998), and CLaMS (McKenna et al., 2002). Exceptions are the LG models ATTILA (Atmospheric Tracer Transport in a LAgrangian model; Reithmeier and Sausen, 2002) and CLaMS, which has recently been coupled to the global chemistry–climate model EMAC (Hoppe et al., 2014). Describing the transport of tracers with an LG transport scheme has advantages compared to an Eulerian transport method: mass conservation (not in CLaMS) and the absence of numerical diffusion. These advantages become most important if tracer distributions are inhomogeneous with strong vertical or horizontal gradients (Stenke and Grewe, 2005), which ought to be smoothed by physical and not by numerical diffusion processes.
ATTILA has already been used to study the advantage of Lagrangian water vapour and cloud water transport on the model climate (Stenke et al., 2008). The results show reduced and thus more realistic water vapour concentrations in the lowermost extratropical stratosphere, a steeper meridional water vapour gradient in the subtropics, and a reduced cold bias around the tropopause near the poles. Furthermore, ATTILA had been used within studies of the climate impact of aviation and climateoptimised air traffic routing (Grewe et al., 2014a, b).
In this study we introduce the extended and improved LG advection scheme ATTILA, which has been parallelised, modularised, and rewritten as a submodel for EMAC (Jöckel et al., 2010). ATTILA was originally developed by Reithmeier and Sausen (2002). We implemented an LG convection scheme and a diabatic vertical velocity formulation, which can be selected instead of the standard kinematic vertical velocity. The need for these two physical improvements is due to the following reasons.
First, the largescale transport of trace species is sensitive to the selected vertical velocity scheme (Eluszkiewicz et al., 2000). Therefore, several studies recommend the use of a diabatic vertical velocity for the representation of LG transport in the tropical tropopause layer and the stratosphere (Eluszkiewicz et al., 2000; Ploeger et al., 2010; Hoppe et al., 2016). Specific transport characteristics like the residence time in the tropical tropopause layer (TTL) and the pathways to the stratospheric overworld are simulated more realistically (Ploeger et al., 2010; Hoppe et al., 2014, 2016) with a diabatic vertical velocity. Typically, kinematic velocities are calculated as a residual from the horizontal flux divergence using the continuity equation. Because horizontal velocities are 2 orders of magnitude larger than the vertical velocity, kinematic velocities show up as rather noisy.
Second, convective transport is an important fast vertical transport process for trace species in the troposphere, and tracer distributions are sensitive to the convection parameterisation (Mahowald et al., 1995; Tost et al., 2006; Erukimova and Bowman, 2006; Zhang et al., 2008). The vertical tracer distribution depends on the accuracy of transport from the boundary layer, where the chemical species are emitted, into the free troposphere. Two LG transport schemes are known to use a convection parameterisation: the CLaMS transport model considers convection by using the moist Brunt–Väisälä frequency parameterisation to include the effects of vertical instability on the related convection (Tao, 2016; Konopka et al., 2018). In the FLEXPART transport model (Stohl et al., 1998) the convection scheme relies on the ECMWF gridscale temperature and humidity and provides a matrix for vertical convective particle displacement (Seibert et al., 2002; Forster et al., 2007).
In the former (nonparallelised) version of ATTILA, convective tracer tendencies were calculated in gridpoint space and then transformed onto the parcels (Reithmeier and Sausen, 2002). This transformation, however, is not mass conserving. Moreover, parcel trajectories do not follow convective updrafts and downdrafts. This is a drawback with respect to the analysis of trajectories, which were subject to convective uplift, and the motivation to incorporate an LG convection scheme in ATTILA. Besides the transport of parcels, mixing of compounds between adjacent parcels is an important process that reduces gradients of trace gases horizontally and vertically. Physically, the character of turbulence in the atmosphere (due to wind shear or buoyancy) controls the degree of mixing. An LG model that successfully uses a physical parameterisation for mixing based on the atmospheric flow deformation is CLaMS (McKenna et al., 2002; Konopka et al., 2004; Riese et al., 2012). However, our parameterisation of mixing, realised in the new submodel LGTMIX, is so far only based on two parameters: one for the troposphere and one for the stratosphere, and it represents local isotropic turbulent mixing. This concept was already successfully applied by Reithmeier and Sausen (2002). However, LGTMIX is written to more easily allow for the incorporation of more physically sound mixing parameterisations in the future.
In Sect. 2 we shortly repeat the main concepts of ATTILA, which were published in detail by Reithmeier and Sausen (2002), and we introduce the application and extensions of the MESSy infrastructure and the concept of the calculation of random numbers in a parallel computing environment. Additionally, we describe the new LG convection parameterisation and the diabatic velocity of the new ATTILA version. The turbulent mixing of compounds between the parcels (submodel LGTMIX), the extended diagnostics (submodel LGVFLUX), the extensions to the submodel DRADON to handle ^{222}Rn in LG representation, and the submodel LGGP, calculating the transformations between GP and LG representation, are also described in Sect. 2. The observational data for comparison are described in Sect. 3. Section 4 describes the model simulations performed with ATTILA coupled to the global chemistry–climate model EMAC. The evaluation of the LG simulations is presented in Sect. 5. We compare the LG simulation results with observations and also with EMAC (GP) simulations, which were already evaluated by Jöckel et al. (2016).
2.1 EMAC – a MESSyfied global chemistry–climate model
The ECHAM/MESSy Atmospheric Chemistry (EMAC) model is a numerical chemistry and climate simulation system that includes submodels describing tropospheric and middle atmosphere processes and their interaction with oceans, land, and human influences (Jöckel et al., 2010). It uses the second version of the Modular Earth Submodel System (MESSy2) to link multiinstitutional computer codes. The core atmospheric model is the fifthgeneration European Centre Hamburg general circulation model (ECHAM5; Roeckner et al., 2006). For the present study we applied EMAC (ECHAM5 version 5.3.02, MESSy version 2.53.0 in the T42L47MAresolution, i.e. with a spherical truncation of T42, corresponding to a quadratic Gaussian grid of approximately 2.8^{∘} by 2.8^{∘} in latitude and longitude) with 47 vertical hybrid pressure levels up to 0.01 hPa (middle atmosphere). The applied model setup comprised the submodels listed in Table 1.
Dietmüller et al. (2016)Reithmeier and Sausen (2002)Roeckner et al. (2006, and references therein)Dietmüller et al. (2016)Tost et al. (2006)Tost (2006)Jöckel et al. (2010)Roeckner et al. (2006, and references therein)Roeckner et al. (2006, and references therein)Sander et al. (2014)Kerkweg et al. (2006)Roeckner et al. (2006, and references therein)Roeckner et al. (2006, and references therein)Jöckel et al. (2008)Giogetta and Bengtson (1999); Jöckel et al. (2006)Dietmüller et al. (2016)Nissen et al. (2007); Dietmüller et al. (2016)Roeckner et al. (2006, and references therein)Kerkweg et al. (2006)Jöckel et al. (2006)Jöckel et al. (2010)The following list gives an overview of the modified and newly developed routines, which are presented in more detail in the following sections.
 a.
Modifications and extensions of physical processes included
 –
additional subroutines for ATTILA to describe Lagrangian convection,
 –
a formulation of vertical movement of air parcels in ATTILA based on the diabatic vertical velocity,
 –
a new submodel (LGTMIX) to calculate the mixing of compounds in Lagrangian representation, and
 –
expansion of the submodel DRADON to account for the emission and decay of ^{222}Rn for the new Lagrangian representation of tracers.
 –
 b.
New diagnostic and infrastructure submodels included
 –
a new submodel for the infrastructure, such as for the calculation of random numbers in a parallel environment,
 –
a subsubmodel that hosts the basic transformation routines needed in ATTILA to convert variables from gridpoint to Lagrangian representation and vice versa (ATTILA_TOOLS),
 –
a new submodel that uses ATTILA_TOOLS to calculate the transformation of userspecified variables between Lagrangian and gridpoint space and vice versa (LGGP) for the output, and
 –
a new submodel to diagnose the vertical fluxes through horizontal surfaces (LGVFLUX).
 –
2.2 Submodel ATTILA: Atmospheric Tracer Transport In a LAgrangian model
ATTILA is a Lagrangian tracer transport scheme, now including LG convection, which can optionally be selected to transport tracers in Lagrangian representation in addition to the standard fluxform semiLagrangian (FFSL) scheme (Lin and Rood, 1996) for tracers in GP representation.
ATTILA runs online as a submodel within EMAC. A former version of ATTILA has been described in detail by Reithmeier and Sausen (2002). The main concepts of ATTILA are shortly repeated in this section (timestepping procedure, interpolation methods, initialisation) and complemented by new infrastructure (random number generator, parallelisation, transformation and transposition methods), as well as new physical (air parcel mixing, Lagrangian convection) and diagnostic submodels.
In ATTILA the atmospheric mass is divided into single mass packets, which have an equal air mass loading but no volume. The parcels are regarded as centroids when they are advected with the wind field provided by the spectral dynamical core of EMAC. The number of parcels within the atmosphere is only limited by the available computational resources. A typical choice is an average of three parcels per EMAC grid box, similar to what was documented by Reithmeier and Sausen (2002). However, the actual number of parcels per grid box may vary between zero and 10, depending on the vertical and horizontal size of the grid box.
2.2.1 Model infrastructure
To enable ATTILA in a distributed memory parallel environment (e.g. applying a messagepassing interface standard, MPI) we chose to follow a domain cloning approach. Whereas the base model EMAC follows a classical horizontal domain decomposition approach for distributed memory parallelisation, we distribute the global number N of ATTILA air parcels, which keep their identity throughout a simulation, (almost) equally among the parallel tasks (index i):
where p is the number of parallel tasks and n_{i} the number of parcels bound to task physical and not by numerical i. Note that n_{i}=n_{j} for all (i,j), except for $i=p\mathrm{1}$, depending on whether N is divisible by p or not.
During the simulation, each parcel keeps being bound to its initial task. Since all parcels on each task move around the entire globe with time, it is necessary to provide the required input variables to drive ATTILA (such as the wind velocity vector from EMAC) as global fields (i.e. by cloning of the global domain of these variables). The subroutines for data transpositions between parallel decomposed grid points and corresponding cloned global variables have been added to the MESSy infrastructure submodel TRANSFORM.
To facilitate the exchange of Lagrangian objects between Lagrangianenabled submodels as socalled channel objects (see Jöckel et al., 2010, for a detailed explanation of the MESSy infrastructure submodel CHANNEL), we define a new representation (see Jöckel et al., 2010, for a detailed explanation) of rank 1, global dimension length N, and local (i.e. task specific) dimension length n_{i}. The corresponding MPIbased gather and scatter routines for serial NetCDF I/O have been added to the MESSy infrastructure submodel TRANSFORM. This new representation, named LG_ATTILA, is used by the Lagrangian submodels to define their specific Lagrangian objects.
For tracers we further define two additional tracer sets (see Jöckel et al. (2008) for a detailed explanation), one (“tracer_lg”) in the new Lagrangian representation to handle the Lagrangian tracers and one “tracer_lggp” in gridpoint representation. The latter is solely used to transform the Lagrangian tracers into gridpoint space for output and further analyses.
Subroutines to transform and transpose variables between Lagrangian representation and (parallel decomposed) gridpoint representation, and vice versa, are collected in a specific toolbox module named ATTILA_TOOLS. This also comprises specific subroutines for the transformation of gridpoint emission fluxes into Lagrangian tracer tendencies. For the latter, four options are implemented: the emitted mass from a grid cell is distributed in one of the following ways.

Evenly among all LG parcels in that grid box. In the case that there is no parcel at a given time in that grid box, the mass is stored and accumulated over time and eventually released into the next parcel(s) passing by.

Evenly among all LG parcels in the lowest grid box of the boundary layer with at least one parcel in it. In the case that the entire column in the boundary layer is empty (i.e. no parcels) at a given time, the mass is stored, accumulated, and eventually released to the next parcel as in (1).

Among all parcels in the boundary layer. However, it is weighted with a linear, negative vertical gradient. The treatment of empty boundary layer columns is as in (1) and (2).

Evenly among all parcels in the boundary layer. The treatment of empty boundary layer columns is as in (1), (2), and (3).
ATTILA requires up to four series of pseudorandom numbers, one for the
boundary layer turbulence parameterisation (Sect. 2.2.3), one for the
convection parameterisation (Sect. 2.2.4), one for an envisaged
(but not yet implemented) additional clear air turbulence parameterisation,
and one for particle displacements parameterising a Monte Carlo diffusion
approach. One additional pseudorandom number series is used for the initial
distribution of the parcels in the model atmosphere. These pseudorandom
number series are provided by the MESSy infrastructure submodel RND. RND
provides uniformly distributed pseudorandom numbers between 0 and 1,
calculated with the standard Fortran90 function,
RANDOM_NUMBER
, the Mersenne Twister algorithm (Matsumoto and Nishimura, 1998), or
the Luxury algorithm (Lüscher, 1994; James, 1994). Based on these, RND can also
provide normally distributed random numbers centred around zero using the
Marsaglia polar
method^{1}. The generation of highquality pseudorandom
number series in a parallel environment is not straightforward. Seeding
independent series on each task implies the high risk that these series
will become correlated. Moreover, the result is decomposition dependent; i.e. it
depends on the number of tasks, which is not desirable. One solution is to
seed one common series on one task and to distribute the resulting
pseudorandom numbers to all other tasks. This implies a load imbalance and
requires additional MPI communication, yet for most pseudorandom
number generators it is the only possibility. However, for the Mersenne
Twister^{2} (among others) Haramoto et al. (2008) found an
efficient “jump ahead” facility, i.e. a method to advance the
pseudorandom number generator state vector by $j=a\times {\mathrm{2}}^{b}$ steps (a
and b integer with a>0 and b>0), without the need to harvest all
j pseudorandom numbers. Jumping ahead by numbers not representable in the
form a×2^{b} can be achieved by additionally harvesting r
pseudorandom numbers such that ${j}^{\prime}=a\times {\mathrm{2}}^{b}+r$. This procedure can
nicely be used for a parallel decomposition independent method for the parallel
generation of pseudorandom number series and has been implemented in the
MESSy infrastructure submodel RND. The same pseudorandom number series is
seeded on all tasks, which then jump ahead and harvest independently, i.e.
without additional communication overhead between the tasks. Each task can
jump ahead directly to the chunk of pseudorandom numbers it needs to
harvest. The only prerequisite for this to work is that the number of
required pseudorandom numbers per (each) task is a priori known to all other
tasks. For instance, if for each ATTILA parcel k pseudorandom numbers are
required (e.g. per model time step), in total k×N pseudorandom
numbers need to be harvested, i.e. k×n_{i} for task i. That means
that task i needs to jump ahead by
steps before it can harvest its own chunk of k×n_{i} pseudorandom numbers.
For the simulations analysed below, we used three uniformly distributed pseudorandom number series, all generated with the Mersenne Twister algorithm: for the boundary layer turbulence scheme, for the convection parameterisation, and for the initial placement of the Lagrangian parcels. The Monte Carlo diffusion was switched off.
2.2.2 Advection
For every time step (in our simulations: Δt=600 s), the parcels are advected by the threedimensional wind field using a fourthorder Runge–Kutta method. The wind field is interpolated on the parcel positions by linear interpolation horizontally (i.e. on the latitude–longitude grid) and by cubic Hermite interpolation vertically. The initialisation of the positions in the atmosphere is carried out randomly so that the number of parcels corresponds to the mass of the respective model layer.
In the vertical direction we may use either ηcoordinate vertical velocities ($\mathit{\eta}=\frac{p}{{p}_{\mathrm{0}}}$, $\dot{\mathit{\eta}}$kinematic velocity) calculated from the horizontal flux divergence using the continuity equation or isentropic coordinates ξ , where the vertical velocities $\dot{\mathit{\xi}}$ are calculated from the EMAC diabatic heating rates (diabatic velocity). The kinematic velocity is provided by default from EMAC, whereas the diabatic velocity was newly implemented similar to Eluszkiewicz et al. (2000) and Hoppe et al. (2016). In our notation,
with θ being the potential temperature and f being defined as
with $\mathit{\kappa}={R}_{\mathrm{v}}/{c}_{\mathrm{p}}p$ and p being the atmospheric pressure; p_{r} is the atmospheric pressure of the climatological tropopause, and c_{p} is the heat capacity at constant pressure. It characterises the transition from a pure θcoordinate system to the ξcoordinate system. The standard surface pressure is p_{0}=1013.25 hPa, and p_{s} is the actual surface pressure.
The vertical velocity in this coordinate system is defined as the time derivative of Eq. (3):
The diabatic vertical velocity in the troposphere $\dot{\mathit{\xi}}$ for p>p_{r} appears as a mixed velocity between pure diabatic $\dot{\mathit{\theta}}$ and kinematic velocity in the troposphere according to Eq. (6). Only in the stratosphere is $\dot{\mathit{\xi}}$ a pure diabatic velocity.
2.2.3 Turbulence
Every parcel located within the planetary boundary layer (PBL) is randomly displaced in the vertical direction within the corresponding grid cell. This stochastic mixing represents the boundary layer convective mixing process. The boundary layer height is calculated outside of ATTILA within the submodel TROPOP.
2.2.4 Convection
The LG convection scheme uses the mass fluxes of the standard gridbox convection scheme in EMAC (submodel CONVECT) to calculate the convective parcel movement. Therefore, we will first shortly introduce the convection scheme of EMAC (Tiedtke, 1989; Nordeng, 1994) because the LG convection scheme is based on it. Convection in the standard convection scheme is initiated when the convergence of moisture in a vertical column of the atmosphere exceeds a certain threshold value and a convectively unstable layer exists. Three types of convection are distinguished: deep convection occurs if moisture convergence through advection and evaporation at the surface takes place. Shallow convection occurs if moisture convergence is only by evaporation at the surface, and midlevel convection occurs if the criteria of deep and shallow convection are not fulfilled but 90 % relative humidity is reached within the planetary boundary layer.
Convection is parameterised by dividing a vertical column into an area of updraft (superscript u), downdraft (superscript d), and an area of compensating motion in the environment (superscript e). Convective transport in EMAC is parameterised only in the vertical direction as a divergence of the tracer mass fluxes F^{u}=M^{u}X^{u}, F^{d}=M^{d}X^{d}, and F^{e}=M^{e}X^{e}:
where $\stackrel{\mathrm{\u203e}}{X}$ is the tracer mass mixing ratio, M is the mass flux, $\stackrel{\mathrm{\u203e}}{\mathit{\rho}}$ is the air density, $\stackrel{\mathrm{\u203e}}{w}$ is the vertical wind component, and z the height. The quantities with an overbar are horizontal averages over the grid box, and the quantities marked with a prime are the horizontal deviations from the respective gridbox mean variables. M^{u}≥0 , M^{d}≤0, and M^{e} are the mass fluxes of air for updraft, downdraft, and the environment, respectively.
The change in mass fluxes with height is dependent on entrainment and detrainment fluxes.
E^{u} (E^{d}) comprises the entrainment (detrainment) rates due to turbulent exchange of mass through cloud edges, and for the updraft E^{u} only, it implies the organised inflow associated with largescale moisture convergence in cases of deep or midlevel convection. Accordingly, the detrainment rates include the turbulent exchange in updraft, downdraft, and, for the updraft only, the organised outflow at cloud top.
The corresponding tracer mass fluxes are as follows.
The calculation of the convective transport of tracer mass starts with the determination of the type of convection (deep, shallow, midlevel). According to the estimated convective available potential energy (CAPE) the mass flux at cloud base is calculated. Further details of the calculation of the mass fluxes are described by Tiedtke (1989) and Nordeng (1994).
In our LG convection scheme air parcels can follow the updraft, downdraft, or the compensating motion in the environment at a grid column with convection within one time step. The forcing used for the Lagrangian convection scheme is provided by the mass fluxes M of the convection scheme of EMAC for updraft and downdraft, respectively. Probabilities P for each level are calculated from the mass fluxes within a vertical column. Each LG parcel is equipped with a (precalculated) random number (see Sect. 2.2.1). For each parcel, ascend (or descend) in an updraft (downdraft) is applied with probability P. The probability for an air parcel to follow the updraft P is equal to the ratio of the mass of the air parcel moving into the updraft to the mass of air at that level.
If $({M}_{k}{M}_{k+\mathrm{1}})>\mathrm{0}$, which means that the mass flux increases with height, then
with
where p is pressure (hPa), A is area (m^{2}), M represents air mass fluxes ($\mathrm{kg}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$), g is gravity acceleration, and Δt is time step length (s).
If $({M}_{k}{M}_{k+\mathrm{1}})<\mathrm{0}$, i.e the mass flux decreases with height, a negative probability is defined to reflect a situation in which a parcel may leave the updraft due to detrainment. The probability is equal to the ratio of the mass leaving the level to the mass entering the same level from below:
The equations of the probability functions are analogous for the downdraft.
The LG convection scheme strictly conserves local mass because for every time step the number of parcels per grid box after convection equals the number before convection (see n= const. in Fig. 1). Every updraft and downdraft forces a compensating largescale motion of parcels. The probability P for subsidence is not estimated from the mass fluxes provided by EMAC. It is calculated for every layer, depending on the number of parcels that need to subside in order to fulfil the mass conservation for every layer.
2.3 Submodel LGGP: transformation between Lagrangian and Eulerian representation
The submodel LGGP (LaGrangian to Grid Point transformations) performs the transformation of variables from Lagrangian representation to gridpoint representation or vice versa. The variables (channel objects) to be transformed are specified by the user in the &CPL namelist of the submodel.
Transformations of a variable from LG to GP use the information of all parcels in the corresponding grid box and calculate

the sum of this variable over all parcels,

the average of the variable over all parcels,

the standard deviation of the variable over all parcels, or

the average of the variable over all parcels in which the variable is >0.
Grid boxes without parcels are either filled with a constant value (defined by the user in the &CPL namelist) or with the value from a selected gridpoint variable (defined as channel object in the &CPL namelist).
The transformation from GP to LG distributes the variable onto all parcels in the respective grid box, either mass conserving (i.e. with equal share) or uniformly (i.e. with the same value of the GP variable). An example &CPL namelist is shown in the Supplement.
2.4 Submodel LGTMIX: mixing of compounds in Lagrangian representation
The submodel LGTMIX (LaGrangian Tracer MIXing) calculates the exchange of tracer mass between Lagrangian parcels. Each Lagrangian parcel is described by a mathematical point. Its tracer mixing ratio represents a mean over the whole parcel. Turbulence in the ambient air leads to the mixing of air of adjacent parcels. In order to avoid parceltoparcel communication, we define a background mixing ratio $\stackrel{\mathrm{\u203e}}{c}$, with which the parcel can communicate. The background is defined by the mean mixing ratio of the individual parcels c_{i} within one grid box of the EMAC grid:
The altered mixing ratio of the respective parcel is then calculated by ${c}_{\mathrm{i}}^{\mathrm{new}}={c}_{\mathrm{i}}+(\stackrel{\mathrm{\u203e}}{c}{c}_{\mathrm{i}})\phantom{\rule{0.25em}{0ex}}d$, with d being a dimensionless mixing parameter within the range [0,1], which controls the magnitude of the exchange.
The user can specify in the LGTMIX &CPL namelist the mixing parameter d individually for vertical model level ranges defined by two external layer definitions (i.e. external channel objects), such as the boundary layer height (from TROPOP), the tropopause (from TROPOP), or any surface provided by VISO (a diagnostic submodel to diagnose vertically layered 2D isosurfaces in 3D scalar fields and to map 3D scalar fields in GP representation on isosurfaces; Jöckel et al., 2010). The value for each of these layers can either be a constant or a function defined by scaling an external gridpoint variable (channel object) in a given range ($[min,max]$, to be specified by the user) to the interval [0,1]. In the latter case, d can also be time dependent. The d in each of these layers can be scaled further for each tracer individually. An example &CPL namelist is shown in the Supplement. For our simulations discussed below, standard values of $d={\mathrm{10}}^{\mathrm{3}}$ for the troposphere and $d=\mathrm{5}\times {\mathrm{10}}^{\mathrm{4}}$ for the stratosphere follow Reithmeier and Sausen (2002). Additional tracers for diagnostic purposes have been simulated without mixing in the stratosphere (i.e. d scaled by 0.0) and with doubled mixing strength (i.e. d scaled by 2.0) in the stratosphere.
2.5 Submodel LGVFLUX: diagnostic of vertical fluxes through horizontal surfaces
The submodel LGVFLUX is a useful tool to calculate online vertical mass fluxes through horizontal surfaces. Mass fluxes through a twodimensional surface (e.g. isentropic surface, potential vorticity isosurface, pressure level) are calculated by analysing the movement of LG particles through these surfaces (upward or downward) and summing over all particles which cross the surface per unit time and area:
with F_{sfc} being the mass flux through the horizontal surface (indicated by the subscript “sfc”), and m_{i} the mass of an LG parcel that is transported through the surface with area A in time Δt (i.e. the model time step length). For air, the mixing ratio c_{i} is 1.0, and for tracer mass fluxes c_{i} is the corresponding tracer mixing ratio. In order to avoid summation over fast, reversible transitions, each surface definition is associated with a minimum residence time each parcel needs after crossing the surface. For taking into account this minimum residence time, each parcel is equipped with a clock to directly measure its transit time. If the parcel crosses a selected horizontal surface, its clock is started and will be reset only if the parcel moves across the surface into the opposite direction. Thus, these “clocks” represent the transit time since passing through a specific surface. In Sect. 5 we present results calculated with this new tool to diagnose the stratosphere–troposphere exchange of air mass and to estimate the age of air (AoA) and the AoA spectra from the transit times in the stratosphere. An example &CPL namelist is shown in the Supplement.
2.6 DRADON
The submodel DRADON (diagnostic radon tracer in GP space (see Sect. 6.1 in Jöckel et al., 2010) has been expanded to handle ^{222}Rn and its decay products as tracers in Lagrangian representation (see Sect. 2.2.1). Here, we simulate ^{222}Rn with a constant ^{222}Rn source of 10 000 $\mathrm{atoms}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$ over icefree land (zero elsewhere) and a decay with a halflife of 3.8 d as the only ^{222}Rn sink.
Further, for the transformation of emission fluxes in GP space into Lagrangian tracer tendencies, the new routines of ATTILA_TOOLS (see Sect. 2.2.1) have been used. As an emission method we selected option 2 (see Sect. 2.2.1); i.e. we put the emitted ^{222}Rn mass into parcels which reside lowest in the boundary layer. An example &CPL namelist is shown in the Supplement.
3.1 ^{222}Rn
^{222}Rn has been frequently used for an evaluation of largescale and convective transport processes (Dentener et al., 1999; Denning et al., 1998; Mahowald et al., 1995, 1997; Jacob et al., 1997; Gupta et al., 2004; Zhang et al., 2008; Jöckel et al., 2010), particularly due to its short lifetime. We selected ^{222}Rn measurements with an annual cycle at different sites over the globe as published by Zhang et al. (2008) and vertical profiles from Kritz et al. (1998) and Zaucker et al. (1996). ^{222}Rn is emitted from land surfaces due to a radioactive decay of radium in soils. ^{222}Rn has a characteristically short radioactive halflife (τ_{½}=3.8 d). For the evaluation of the simulated ^{222}Rn distribution we use monthly mean surface values from 18 stations worldwide and selected vertical profiles. The large set of monthly mean surface values was collected from the literature by Zhang et al. (2008) and is used here for comparison. Observations of the vertical distribution of ^{222}Rn are rare, especially if they cover more than the boundary layer. We use two different data sets of vertical profiles for comparison, as outlined below.
Kritz et al. (1998) used flights from the Kuiper Airborne Observatory in summer 1994 to achieve a representative selection of ^{222}Rn measurements in the free troposphere. The flights were made from 3 June until 16 August 1994 around Moffett Field in California (37.4^{∘} N, 122^{∘} W), where 11 single profiles could be realised with a vertical resolution of 1 km.
Zaucker et al. (1996) compiled a data set from nine flights in August 1993 from cities in Nova Scotia and New Brunswick on the east coast of Canada to the western North Atlantic Ocean during the North Atlantic Regional Experiment (NARE) Intensive. The vertical height of the measurements is restricted from the surface to about 5.5 km.
3.2 Age of air
Mean age of air (AoA) is a common metric to quantify the overall capabilities of a global model to simulate stratospheric transport. It describes the transit time of air parcels in the stratosphere (Hall and Plumb, 1994). AoA is calculated (in the model and observations) from an inert tracer with linearly increasing boundary conditions at the surface. AoA at a certain grid point in the stratosphere is then calculated as the time lag between the local tracer mixing ratio and the mixing ratio at a reference point (e.g. the boundary layer in the tropics). Inert tracers from observations are anthropogenically emitted sulfur hexafluoride (SF_{6}) (Stiller et al., 2012; Haenel et al., 2015) and CO_{2} (Engel et al., 2009; Andrews et al., 2001). Both will be used for comparison.
We performed two identical simulations with EMAC–ATTILA with respect to the climate: one uses the kinematic vertical velocity to drive the Lagrangian parcels, and the other uses the diabatic vertical velocity. The horizontal velocity remains equal in both simulations. EMAC was operated in T42L47MA resolution with 47 levels up to 0.01 hPa (see Sect. 2). We simulated the years 1950 to 2010 with prescribed sea surface temperatures (SSTs) from the global data set HadISST (available from http://www.metoffice.gov.uk/hadobs/hadisst/, last access: 6 May 2019; Rayner et al., 2003) similar as for the RC1base08 freerunning simulation (see ESCiMo project description by Jöckel et al., 2016). In contrast to the simulations within the ESCiMo project, our simulations do not simulate interactive chemistry; however, monthly averages of radiatively active substances (CO_{2}, O_{3}, N_{2}O, CF_{2}Cl_{2}, CFCl_{3}) have been prescribed from RC1base08. Methane was initialised as in RC1base08 and the same pseudoemission time series (by Newtonian relaxation at the surface with the submodel TNUDGE) has been applied. However, methane oxidation and its contribution to stratospheric water vapour were treated in a simplified manner (with the submodel CH4): the oxidation educts OH, Cl, and O^{1}D have been prescribed as monthly averages from RC1base08. The photolysis rate J(CH_{4}) was calculated with the submodel JVAL.
ATTILA was initialised with 1.15×10^{6} parcels, which in sum represent the total mass of the atmosphere. The parcels were initially positioned according to the mass distribution in grid space. The results of the two simulations are further denoted as

GP for the results of the gridpoint simulation (EMAC; note that these are identical in both simulations),

LG(diab) for the results of EMAC–ATTILA with diabatic vertical velocity, and

LG(kin) for the results of EMAC–ATTILA with kinematic vertical velocity.
The LG parcels are equipped with tracers with different properties.

^{222}Rn is a commonly used tracer to study vertical transport into the upper troposphere due to its characteristically short radioactive halflife of (τ_{½}=3.8 d). In our simulations it is emitted at the surface with an emission rate of 10 000 $\mathrm{atoms}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}{\mathrm{s}}^{\mathrm{1}}$ over icefree land surfaces between 90^{∘} N and 90^{∘} S.

SF_{6}_AoA and SF_{6}_AoAc are inert synthetic tracers. They differ with respect to the surface source. Both are nudged by Newtonian relaxation at the lowest model layer towards a linearly increasing mixing ratio. Note that for SF_{6}_AoA the linearintime increasing mixing ratio is latitude dependent. Using a spatially constant surface source (SF_{6}_AoAc) has the advantage that concentration differences in the atmosphere cannot have their origin in the distribution of surface sources.

SF_{6}_AoA_nm has the same properties as the tracer SF_{6}_AoA; however, the interparcel mixing was set to zero (see Sect. 2.4). Hence, it can be used in comparison to SF_{6}_AoA to study the influence of local mixing between adjacent parcels on the global AoA distribution.
In the previous sections, we described a comprehensively updated version of the LG tracer transport scheme ATTILA, including a new LG convection scheme and the option to use a diabatic instead of the standard kinematic vertical velocity. In this section, we evaluate ATTILA by comparing the simulated ^{222}Rn and SF_{6}_AoA tracer distributions with observations and with EMAC results, i.e. from the GP space.
5.1 Simulation of ^{222}Rn
Jöckel et al. (2010) showed that the EMAC model (version 2.40) is able to realistically simulate the ^{222}Rn distribution. We therefore assume here that our GP simulation with the EMAC model (now version 2.53) simulates ^{222}Rn similarly. A comparison with observations will follow in the next section. The intercomparison of our new simulations of ^{222}Rn between GP and LG space (Fig. 2) shows the largest values of ^{222}Rn in the northern hemispheric boundary layer (the lowest three model layers). The large maxima north of 10^{∘} N and around 30^{∘} S are related to surface emissions from the large continents (Fig. 3). The small local maximum at 80^{∘} S is related to the surface edge of the Antarctic continent, where small land areas in the land–sea mask generate local ^{222}Rn emissions. In the boundary layer (north of 40^{∘} N) the zonal mean ^{222}Rn values for LG(diab) are smaller than for LG(kin) and GP (Fig. 2). In contrast, south of 40^{∘} N the LG(diab) results are closer to the GP simulation. Part of the difference between the GP and the LG simulation is the uptake of emissions, which depends on the number of LG parcels present in the lowest model levels (see also Sect. 2.6) and their relative horizontal distribution (since the source is only over land). As expected, over the oceans (the remote regions) the ^{222}Rn values are relatively small. At 100 hPa in Fig. 2 the largest ^{222}Rn values occur in the tropics. Here, GP and LG(diab) simulation results are in close agreement. LG(kin) simulates smaller values. The differences between LG(diab) and LG(kin) are related to the different vertical velocity scheme because this is the only difference between the two LG schemes. In the higher latitudes the 100 hPa level is already in the stratosphere and ^{222}Rn has largely decayed due to its short halflife and the relatively long transport times in the stratosphere.
5.1.1 Annual cycle of ^{222}Rn at the surface layer
We use ^{222}Rn in situ measurements at the surface layer of 18 stations distributed worldwide as they were published by Zhang et al. (2008). The model results are longterm averages of ^{222}Rn in the surface layer over the years 1960–2000. They are horizontally linearly interpolated to the respective location of the observations. Six stations (Crozet, Bermuda, Amsterdam Island, Kerguelen, Dumont, and Mauna Loa) are far away from the continents and show the effect of longrange transport (^{222}Rn lower than 1000 mBq m^{−3}; STP). A further six stations (Socorro, Cincinnati, Para, Puy de Dôme, Beijing, and Hohenpeißenberg) are located on continents in the vicinity of the ^{222}Rn sources. Finally, six stations (Gosan, Hong Kong, Cape Grim, Livermore, Bombay, and Mace Head) are located at coastal sites influenced by the prevailing wind direction either from the sea or the continent. A detailed comparison between model results (horizontal axis) and measurements (vertical axis) is shown in Fig. 4. A total of 12 monthly values of the mean annual cycle of GP, LG(diab), and LG(kin) are presented as data points with different colours for each station. The results are similar to those of Jöckel et al. (2010, their Fig. 14). This cannot be taken for granted because Jöckel et al. (2010) used a nudged simulation and a model setup with 90 vertical levels, whereas in our study the simulation is free running and our model setup has 47 vertical levels. Two stations are for all 12 monthly mean values out of the thick dashed area in Fig. 4: Beijing and Kerguelen (southern Indian Ocean). Beijing's ^{222}Rn values are too low and the Kerguelen values are too high as simulated with all models. In Jöckel et al. (2010) the measured ^{222}Rn values for Kerguelen are not captured in the simulation. The simulated ^{222}Rn values are strongly dependent on the local wind direction at the surface of the measurement sites, especially for the coast and the remote regions over sea.
5.1.2 Vertical profiles of ^{222}Rn
The NARE campaign took place in the vicinity of Nova Scotia and Brunswick, Canada, in August 1993. Data were sampled over the ocean and over the continent. We used the simulation data of a climatological mean August (1960–2000) averaged over the region where the flights took place (60–70^{∘}W and 41–46^{∘}N). The triangles in Fig. 5 show the measured ^{222}Rn in the atmosphere for several flights. The thick dashed curve is a spline interpolation on the grid of all flights. The ^{222}Rn mixing ratio decreases with height up to about 3 km and remains around 10^{−24} mol mol^{−1}. Both LG simulation results agree with the observations and are close to the GP results. The measurements of the Moffett campaign (Fig. 6) show a relatively large scatter of 11 single profiles in the free troposphere. The simulated ^{222}Rn profiles in that region were selected as a climatological mean for the months June to August (1960–2000) and capture the observations quite well. Observed ^{222}Rn emissions are highly variable due to the dependence on the physical characteristics of the soil (Gupta et al., 2004; Karstens et al., 2015). Therefore, a certain spread between observations and model results is expected and acceptable.
5.2 Age of air
The calculation of AoA is performed in two different ways: by a socalled clock tracer (a linearintime increasing tracer like SF_{6}) or directly by a clock on a parcel (LG clock; Sect. 2.5). From a clock tracer, AoA is calculated indirectly by comparing local tracer mixing ratios with reference mixing ratios, e.g. the surface mixing ratio in the tropics (see Sect. 3.2). This concept is applied in the next section for the calculation of mean AoA. However, age spectra are calculated directly from transit times provided by the parcel clocks (e.g. during the transit in the stratosphere). Conceptually, the calculated AoA differs if it was calculated by a clock tracer or by clocks. A clock tracer is subject to interparcel mixing, which is not the case for parcel clocks. Therefore, the mean AoA calculated by the parcel clocks is older than for a clock tracer. However, AoA from a simulated clocktracer distribution can be directly compared to AoA from an observed tracer distribution.
5.2.1 Mean age of air
Mean AoA in the stratosphere is calculated from the SF_{6}_AoA tracer. The transit time is estimated by comparing the tracer mixing ratio in the stratosphere with the SF_{6}_AoA mixing ratio in the tropical boundary layer. Figure 7 shows a comparison of SF_{6}_AoA at 20 km of height between the GP, LG(kin), and LG(diab) results along with AoA derived from satellite observations from MIPAS (Stiller et al., 2012) and in situ measurements of Waugh and Hall (2002) from the SF_{6} tracer, as well as from CO_{2} in situ data (Andrews et al., 2001). GP and LG(diab) show realistic distributions of AoA, although slightly lower AoA than the in situ measurements. MIPAS data are known to overestimate AoA in the polar regions (Stiller et al., 2012). This is attributed to a known sink of SF_{6} in the upper stratosphere that is not accounted for in our simulations. Noticeably, the LG(diab) simulation is closer to the observations and LG(kin) shows up with an age that is too low. For the analysis, we therefore restrict our further evaluation to LG(diab). The mean age distribution (19602010) of LG(diab) is shown in Fig. 8. It confirms the wellknown characteristics of the stratospheric Brewer–Dobson circulation with younger air in the tropical pipe and older air over the poles. Furthermore, the simulated AoA is slightly older in the whole stratosphere compared to GP (Fig. 9), which is most pronounced near the poles below 50 hPa. This difference is attributed to the Eulerian vertical velocity used in the fluxform semiLagrangian transport scheme for the GP simulations, which shows upwelling and downwelling at different high latitudes that are not related to the net tracer transport (Hoppe et al., 2016).
5.2.2 Typical age spectra
AoA spectra are calculated directly from the clock transit times. These LG clocks represent the actual time a parcel resides in the stratosphere after it has crossed the tropopause level. However, these clocks do not “mix their time” with other parcels. Therefore, the resulting spectrum might differ from age spectra calculated from socalled AoA “clock tracers” (Garny et al., 2014; Ploeger et al., 2014; Schoeberl et al., 2005). Typical age spectra for the tropics (20^{∘} N–20^{∘} S) and the poles (70–90^{∘} N, 70–90^{∘} S) between 50 and 0.1 hPa are shown in Fig. 10. The frequency distribution of the LG(diab) simulation is calculated for every month for the years 1990–2010, binned into 0.5year bins, and normalised. We find the largest frequency at a parcel age between 0 and 0.5 years in the tropics and an exponential shape of the spectrum. In contrast, the modal age is roughly 3 years at the poles. The shape of the spectrum is Gaussian with a positive skewness. The seasonal age spectra (selected between 400 and 500 K at 50–70^{∘} N; see Fig. 11) show characteristic multiple local maxima along the time axis. The width of these maxima increases from spring (MAM) and summer (JJA) over autumn (SON) to winter (DJF). The distance of the maxima on the transit time axis is around 1 year. These maxima reflect the different contributions of air masses from the tropics and from high latitudes in the seasonal cycle (Ploeger and Birner, 2016). The seasonal age spectra look qualitatively similar as in Fig. 5 of Ploeger and Birner (2016), although our modal values are about 0.3 years younger. The different modal values between Ploeger and Birner (2016) and our Fig. 11 are probably a consequence of the utilised concept in calculating the AoA. In LG(diab) we use our LG clocks to calculate the transit time in the stratosphere directly if the parcels cross the tropopause level (see Sect. 2.5). Ploeger and Birner (2016) calculated their seasonal spectrum from an AoA clock tracer and the Green's function (Waugh and Hall, 2002) and relate their stratospheric mixing ratios to the tracer mixing ratio in the boundary layer to calculate the AoA. Therefore, the larger modal values in Ploeger and Birner (2016) refer to an additional transit time up to the tropopause level. The effect of these two different concepts on the calculation of AoA in the lower stratosphere is discussed further in Sect. 5.2.3.
5.2.3 Sensitivity of age of air to interparcel mixing
AoA is influenced by the amount of mixing between adjacent parcels. Interparcel mixing can be regarded as a diffusion process leading to a reduction of local AoA gradients. The effect of interparcel mixing makes stratospheric air generally younger (Fig. 12). Stratospheric AoA without interparcel mixing as represented in LG(diab), described as LG(diabnm), is mostly up to 2.5 months older compared to the AoA with standard mixing (see Sect. 2.4) but slightly younger in the tropical lower stratosphere, where mixing with upper tropospheric air becomes important. This “younger air through interparcel mixing” should not be confused with the “ageing by mixing” concept of Garny et al. (2014), Ploeger et al. (2014), Konopka et al. (2015), and Dietmüller et al. (2017). They calculate AoA from the tracer budget equation of AoA and distinguish between the different terms: tendencies due to the residual stratospheric circulation and the tendencies of AoA due to eddy mixing and due to turbulent diffusion. Their concept allows for the separation of the contribution of mixing on the local AoA budget. In contrast, in our study we simply compare a simulated tracer with mixing to a tracer distribution without mixing (perturbation concept). Interparcel mixing in the troposphere (Fig. 13) has only a small and statistically insignificant effect on the simulated AoA because the troposphere is a wellmixed region where parcels often have contact with the surface source. The surface is the only region (in the vertical) where the tracer mixing ratio for the tracer without interparcel mixing can be changed.
5.3 Stratosphere–troposphere exchange
Stratosphere–troposphere exchange (STE) is characterised by a globalscale meridional circulation in which mass is transported upward in the tropics and downward in the extratropics (Holton et al., 1996). We use the new diagnostic submodel LGVFLUX to directly calculate the simulated mass flux through the tropopause in the LG simulation. The LG(diab) simulation captures these typical features (see Fig. 14 as an example for the mass flux through the 380 K isentropic surface) with a net upward flux in the tropics between 30^{∘} S and 30^{∘} N and a net downward flux from 40 to 90^{∘} N and S. Between 30 and 40^{∘} N and S the zonal mean net flux is near zero. Figure 15 shows the annual cycle of the net downward flux through the extratropical 380 K isentropic surface with a maximum in boreal winter for 30–90^{∘} N and boreal summer for 30–90^{∘} S. The annual amplitude between summer and winter is about 6–7 ×10^{9} kg s^{−1} and falls in the range given by Appenzeller et al. (1996) of 6–7 ×10^{9} kg s^{−1}.
5.4 Lagrangian convection statistics
The LG convective parcel movement depends on the calculated mass flux profile (from convection). We analysed the movement of parcels during deep convective events for the year 1997. The analysis of movement shows that within the updraft the largest number of parcels leave the boundary layer and are detrained into the free troposphere up to the tropopause (Fig. 16). The maximum levels of detrainment are between 500 and 600 hPa. Only a few parcels start to follow the updraft above the boundary layer (start level less than 900 hPa). Parcels in the downdraft (Fig. 17) start between 300 and 750 hPa and most parcels are released into the boundary layer (lowest three model layers). Interestingly, three height regions seem to preferably be the starting points for the downdraft: around 400 hPa, between 600 and 700 hPa, and at 850 hPa. The compensating motion in the environment is a movement over a small distance only. The starting levels for subsidence comprise nearly the whole troposphere (Fig. 18). The most frequent movement is one level, but a few parcels are subject to a further downward movement (right side of the diagonal line). However, because the subsidence of parcels depends on a local probability (see Sect. 2.2.4), it is possible that even more parcels subside than originally should. This is then compensated for in the next iteration by the rise of a parcel (left side of the diagonal line). This upward movement of parcels adds a certain amount of unphysical diffusion to the convection that unfortunately cannot be avoided in this model setup.
In this study we described and evaluated the updated LG tracer transport scheme ATTILA. ATTILA was extended with an LG convection scheme and a formulation of diabatic vertical velocity. We implemented a submodel to describe interparcel mixing, which has so far been set up with one parameter for the troposphere and one for the stratosphere. Moreover, the new submodel allows us to easily implement more physically sound mixing parameterisations. New infrastructure submodels which simplify the transformation between GP and LG space, the provision of random numbers in a parallel environment, and diagnostic submodels were developed. We performed two simulations from 1950 to 2010, both resulting in the same meteorological sequence in GP. The simulations differ only with respect to the vertical velocity used for the LG model: one with a diabatic LG(diab) and one with the standard kinematic vertical velocity LG(kin). The annual cycle of the two LG simulations of ^{222}Rn in the surface layer is in accordance with observations of a large number of stations and of comparable quality with a former nudged simulation in grid space. Vertical profiles of ^{222}Rn measured during the NARE and Moffett campaigns agree within an acceptable spread with our simulations. We expected the largest improvement of our results with respect to the simulation of AoA in the stratosphere in the LG(diab) simulation. Indeed, AoA in LG(diab) shows the best agreement with observations. Moreover, AoA spectra and the troposphere–stratosphere exchange are realistically simulated in LG(diab).
In a next step, we plan to parameterise phase changes of water vapour due to convection and cloud development on the Lagrangian parcels. Then, ATTILA will allow us to study convective and largescale water vapour transport consistently in convective regions and to assess the convective contribution to the stratospheric water vapour budget.
The Modular Earth Submodel System (MESSy) is continuously further developed and applied by a consortium of institutions. The usage of MESSy and access to the source code is licenced to all affiliates of institutions which are members of the MESSy Consortium. Institutions can become a member of the MESSy Consortium by signing the MESSy Memorandum of Understanding. More information can be found on the MESSy Consortium Website (http://www.messyinterface.org, last access: 6 May 2019). The code presented here is based on MESSy version 2.53.0 and will be available in the next official release (version 2.55.0). The data from the simulations will be provided by the authors on request.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1219912019supplement.
SB and PJ did the implementation, performed the simulations, analysed the results, and wrote the paper.
The authors declare that they have no conflict of interest.
The data on the annual cycle of ^{222}Rn at several locations in the
surface layer were kindly provided by Kai Zhang. The data on the vertical
radon profiles are kindly provided by Holger Tost. The simulations were
performed at the LeibnizRechenzentrum in Garching, Germany. The data were
analysed with the interactive computer visualisation and analysis environment
ferret. This work was funded by the Deutsche Forschungsgemeinschaft (DFG)
within the projects LAWA (Lagrangesche Simulation des globalen
atmospärischen Wasserkreislaufs) under the grant EG 40/241, the DFG
Forschergruppe SHARP (Stratospheric Change and its Role for Climate
Prediction), the DLRProject KliSAW (Klimarelevanz von atmosphärischen
Spurengasen, Aerosolen und Wolken), and the HelmholtzGemeinschaft e.V. (HGF)
Project
“Advanced Earth System Modelling Capacity (ESM)”.
We thank Oliver Reitebuch for the internal review of an early version of the
paper.
The article processing charges for
this openaccess
publication were covered by a Research
Centre of the Helmholtz Association.
This paper was edited by Havala Pye and reviewed by two anonymous referees.
Andrews, A. E., Boering, K. A., Daube, B. C., Wofsy, S. C., Loewenstein, M., Jost, H., Podolske, J. R., Webster, C. R., Herman, R. L., Scott, D. C., Flesch, G. J., Moyer, E. J., Elkins, J. W., Dutton, G. S., Hurst, D. F., Moore, F. L., Ray, E. A., Romashkin, P. A., and Strahan, S. E.: Mean ages of stratospheric air derived from in situ observations of CO_{2}, CH_{4}, and N_{2}O, J. Geophys. Res.Atmos., 106, 32295–32314, https://doi.org/10.1029/2001JD000465, 2001. a, b, c
Appenzeller, C., Holton, J. R., and Rosenlof, K. R.: Seasonal variation of mass transport across the tropopause, J. Geophys. Res., 101, 15071–15078, https://doi.org/10.1029/96JD00821, 1996. a
Collins, W. J., Derwent, R. G., Johnson, C. E., and Stevenson, D. S.: A comparison of two schemes for the convective transport of chemical species in a Lagrangian global chemistry model. Q. J. Roy. Meteorol. Soc., 128, 9911009, 2002. a
Denning, A. S., Holzer, M., Gurney, K. R., Heimann, M., Law, R. M., Rayner, P. J., Fung, I. Y., Fan, S.M., Taguchi, S., Friedlingstein, P., Balkanski, Y., Taylor, J., Maiss, M., and Levin, I.: Threedimensional transport and concentration of SF6 A model intercomparison study (TransCom 2), Tellus B, 51, 266–297, 1998. a
Dentener, F., Feichter, J., and Jeuken, A.: Simulation of the transport of Rn using online and offline global models at different horizontal resolutions: A detailed comparison with measurements, Tellus B, 51, 573–602, https://doi.org/10.3402/tellusb.v51i3.16440, 1999. a
Dietmüller, S., Jöckel, P., Tost, H., Kunze, M., Gellhorn, C., Brinkop, S., Frömming, C., Ponater, M., Steil, B., Lauer, A., and Hendricks, J.: A new radiation infrastructure for the Modular Earth Submodel System (MESSy, based on version 2.51), Geosci. Model Dev., 9, 2209–2222, https://doi.org/10.5194/gmd922092016, 2016. a, b, c, d
Dietmüller, S., Garny, H., Plöger, F., Jöckel, P., and Cai, D.: Effects of mixing on resolved and unresolved scales on stratospheric age of air, Atmos. Chem. Phys., 17, 7703–7719, https://doi.org/10.5194/acp1777032017, 2017. a
Engel, A., Möbius, T., Bönisch, H., Schmidt, U., Heinz, R., Levin, I., Atlas, E., Aoki, S., Nakazawa, T., Sugawara, S., Moore, F., Hurst, D., Elkins, J., Schauffler, S., Andrews, A., and Boering, A.: Age of stratospheric air unchanged within uncertainties over the past 30 years, Nat. Geosci., 2, 28–31, 2009. a
Eluszkiewicz, J., Hemler, R. S., Mahlman, J. D., Bruhwiler, L., and Takacs, L. L.: Sensitivity of AgeofAir Calculations to the Choice of Advection Scheme, J. Atmos. Sci., 57, 3185–3201, https://doi.org/10.1175/15200469(2000)057<3185:SOAOAC>2.0.CO;2, 2000. a, b, c
Erukhimova, T. and Bowman, K. P.: Role of convection in globalscale transport in the troposphere, J. Geophys. Res., 111, D03105, https://doi.org/10.1029/2005JD006006, 2006. a
Forster, C., Stohl, A., and Seibert, P.: Parameterization of convective transport in a Lagrangian particle dispersion model and its evaluation, J. Appl. Meteorol. Clim., 46, 403–422, 2007. a
Garny, H., Birner, T., Bönisch, H., and Bunzel, F.: The effects of mixing on age of air, J. Geophys. Res.Atmos., 119, 7015–7034, https://doi.org/10.1002/2013JD021417, 2014. a, b
Giorgetta, M. A. and Bengtsson, L.: The potential role of the quasibiennial oscillation in the stratospheretroposphere exchange as found in water vapour in general circulation model experiments, J. Geophys. Res., 104, 6003–6019, 1999. a
Grewe, V., Frömming, C., Matthes, S., Brinkop, S., Ponater, M., Dietmüller, S., Jöckel, P., Garny, H., Tsati, E., Dahlmann, K., Søvde, O. A., Fuglestvedt, J., Berntsen, T. K., Shine, K. P., Irvine, E. A., Champougny, T., and Hullah, P.: Aircraft routing with minimal climate impact: the REACT4C climate cost function modelling approach (V1.0), Geosci. Model Dev., 7, 175–201, https://doi.org/10.5194/gmd71752014, 2014a. a
Grewe V., Champougny, T., Matthes, S., Frömming, C., Brinkop, S., Søvde, A. O., Irvine E. A., and Halscheidt, L.: Reduction of the air traffic's contribution to climate change: A REACT4C case study, Atmos. Environ., 94 616–625, 2014. a
Gupta, M. L., Douglass, A. R., Kawa, S. R., and Pawson, S.: Use of radon for evaluation of atmospheric transport models: sensitivity to emissions, Tellus B, 56, 404–412, https://doi.org/10.1111/j.16000889.2004.00124.x, 2004. a, b
Haenel, F. J., Stiller, G. P., von Clarmann, T., Funke, B., Eckert, E., Glatthor, N., Grabowski, U., Kellmann, S., Kiefer, M., Linden, A., and Reddmann, T.: Reassessment of MIPAS age of air trends and variability, Atmos. Chem. Phys., 15, 13161–13176, https://doi.org/10.5194/acp15131612015, 2015. a
Hall, T. M. and Plumb, R. A.: Age as a diagnostic of stratospheric transport, J. Geophys. Res.Atmos., 99, 1059–1070, 1994. a
Haramoto, H., Matsumoto, M., Nishimura, T., Panneton, F., and L'Ecuyer, P.: Efficient Jump Ahead for F2Linear Random Number Generators, INFORMS J. Comput., 20, 385–390, https://doi.org/10.1287/ijoc.1070.0251, 2008. a
Holton, J. R., Haynes, P. H., McIntyre, M. E., Douglass, A. R., Rood, R. B., and Pfister, L.: Stratospheretroposphere exchange, Rev. Geophys., 33, 403–439, https://doi.org/10.1029/95RG02097, 1995. a
Hoppe, C. M., Hoffmann, L., Konopka, P., Grooß, J.U., Ploeger, F., Günther, G., Jöckel, P., and Müller, R.: The implementation of the CLaMS Lagrangian transport core into the chemistry climate model EMAC 2.40.1: application on age of air and transport of longlived trace species, Geosci. Model Dev., 7, 2639–2651, https://doi.org/10.5194/gmd726392014, 2014. a, b
Hoppe, C. M., Ploeger, F., Konopka, P., and Müller, R.: Kinematic and diabatic vertical velocity climatologies from a chemistry climate model, Atmos. Chem. Phys., 16, 6223–6239, https://doi.org/10.5194/acp1662232016, 2016. a, b, c, d
Jacob, D. J., Prather, M. J., Rasch, P. J., Shia, R.L., Balkanski, Y. J., Beagley, S. R., Bergmann, D. J., Blackshear, W. T., Brown, M., Chiba, M., Chipperfield, M. P., de Grandpré, J., Dignon, J. E., Feichter, J., Genthon, C., Grose, W. L., Kasibhatla, P. S., Köhler, I., Kritz, M. A., Law, K., Penner, J. E., Ramonet, M., Reeves, C. E., Rotman, D. A., Stockwell, D. Z., Van Velthoven, P. F. J., Verver, G., Wild, O., Yang, H., Zimmermann, P.: Evaluation and intercomparison of global atmospheric transport models using ^{222}Rn and other shortlived tracers, J. Geophys. Res., 102, 5953–5970, https://doi.org/10.1029/96JD02955 1997. a
James, F.: RANLUX: A Fortran implementation of the highquality pseudorandom number generator of Lüscher, Comput. Phys. Commun., 79, 111–114, https://doi.org/10.1016/00104655(94)90233X, 1994. a
Jöckel, P., Sander, R., Kerkweg, A., Tost, H., and Lelieveld, J.: Technical Note: The Modular Earth Submodel System (MESSy) – a new approach towards Earth System Modeling, Atmos. Chem. Phys., 5, 433–444, https://doi.org/10.5194/acp54332005, 2005.
Jöckel, P., Tost, H., Pozzer, A., Brühl, C., Buchholz, J., Ganzeveld, L., Hoor, P., Kerkweg, A., Lawrence, M. G., Sander, R., Steil, B., Stiller, G., Tanarhte, M., Taraborrelli, D., van Aardenne, J., and Lelieveld, J.: The atmospheric chemistry general circulation model ECHAM5/MESSy1: consistent simulation of ozone from the surface to the mesosphere, Atmos. Chem. Phys., 6, 5067–5104, https://doi.org/10.5194/acp650672006, 2006. a, b
Jöckel, P., Kerkweg, A., BuchholzDietsch, J., Tost, H., Sander, R., and Pozzer, A.: Technical Note: Coupling of chemical processes with the Modular Earth Submodel System (MESSy) submodel TRACER, Atmos. Chem. Phys., 8, 1677–1687, https://doi.org/10.5194/acp816772008, 2008. a, b
Jöckel, P., Kerkweg, A., Pozzer, A., Sander, R., Tost, H., Riede, H., Baumgaertner, A., Gromov, S., and Kern, B.: Development cycle 2 of the Modular Earth Submodel System (MESSy2), Geosci. Model Dev., 3, 717–752, https://doi.org/10.5194/gmd37172010, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m
Jöckel, P., Tost, H., Pozzer, A., Kunze, M., Kirner, O., Brenninkmeijer, C. A. M., Brinkop, S., Cai, D. S., Dyroff, C., Eckstein, J., Frank, F., Garny, H., Gottschaldt, K.D., Graf, P., Grewe, V., Kerkweg, A., Kern, B., Matthes, S., Mertens, M., Meul, S., Neumaier, M., Nützel, M., OberländerHayn, S., Ruhnke, R., Runde, T., Sander, R., Scharffe, D., and Zahn, A.: Earth System Chemistry integrated Modelling (ESCiMo) with the Modular Earth Submodel System (MESSy) version 2.51, Geosci. Model Dev., 9, 1153–1200, https://doi.org/10.5194/gmd911532016, 2016. a, b
Karstens, U., Schwingshackl, C., Schmithüsen, D., and Levin, I.: A processbased ^{222}radon flux map for Europe and its comparison to longterm observations, Atmos. Chem. Phys., 15, 12845–12865, https://doi.org/10.5194/acp15128452015, 2015. a
Kerkweg, A., Sander, R., Tost, H., and Jöckel, P.: Technical note: Implementation of prescribed (OFFLEM), calculated (ONLEM), and pseudoemissions (TNUDGE) of chemical species in the Modular Earth Submodel System (MESSy), Atmos. Chem. Phys., 6, 3603–3609, https://doi.org/10.5194/acp636032006, 2006. a, b
Konopka, P., Steinhorst, H.M., Grooß, J.U., Günther, G., Müller, R., Elkins, J. W., Jost, H.J., Richard, E., Schmidt, U., Toon, G., and McKenna, D. S.: Mixing and ozone loss in the 19992000 Arctic vortex: Simulations with the 3dimensional Chemical Lagrangian Model of the Stratosphere (CLaMS), J. Geophys. Res., 109, D02315, https://doi.org/10.1029/2003JD003792, 2004. a
Konopka, P., Ploeger, F., Tao, M., Birner, T., and Riese, M.: Hemispheric asymmetries and seasonality of mean age of air in the lower stratosphere: Deep versus shallow branch of the BrewerDobson circulation, J. Geophys. Res.Atmos., 120, 2053–2066, https://doi.org/10.1002/2014JD022429. 2015. a
Konopka, P., Tao, M., Ploeger, F., Diallo, M., and Riese, M.: Tropospheric mixing and parametrization of unresolved convection as implemented into the Chemical Lagrangian Model of the Stratosphere (CLaMS), Geosci. Model Dev. Discuss., https://doi.org/10.5194/gmd2018165, in review, 2018. a
Kritz, M. A., Rosner S. W., and Stockwell, D. Z.: Validation of an offline threedimensional chemical transport model using observed radon profiles, J. Geophys. Res., 103, 8425–8432, 1998. a, b
Lin, S.J. and Rood, R.: Multidimensional fluxform semi Lagrangian transport schemes, Mon. Weather Rev., 124, 2046–2070, 1996. a
Lüscher, M.: A portable highquality random number generator for lattice field theory simulations, Comput. Phys. Commun., 79, 100–110, https://doi.org/10.1016/00104655(94)902321, 1994. a
Mahowald, N. M., Rasch, P. J., and Prinn, R. G.: Cumulus parameterizations in chemical transport models, J. Geophys. Res.Atmos., 100, 26173–26189, 1995. a, b
Matsumoto, M. and Nishimura, T.: Mersenne twister: a 623dimensionally equidistributed uniform pseudorandom number generator, ACM Trans. Model. Comput. Simul., 8, 3–30, 1998. a
Mahowald, N. M., Rasch, P. J., Eaton, B. E., Whittlestone, S., Prinn, R. G.: Transport of ^{222}radon to the remote troposphere using the model of atmospheric transport and chemistry and assimilated winds from ECMWF and the National Center for Environmental Prediction NCAR, J. Geophys. Res.Atmos., 102, 28139–28151, 1997. a
McKenna, D. S., Konopka, P., Grooß, J.U., Günther, G., Müller, R., Spang, R., Offermann, D., and Orsolini, Y.: A new Chemical Lagrangian Model of the Stratosphere (CLaMS): 1. Formulation of advection and mixing, J. Geophys. Res., 107, 4309, https://doi.org/10.1029/2000JD000114, 2002. a, b
Nissen, K. M., Matthes, K., Langematz, U., and Mayer, B.: Towards a better representation of the solar cycle in general circulation models, Atmos. Chem. Phys., 7, 5391–5400, https://doi.org/10.5194/acp753912007, 2007. a
Nordeng, T. E.: Extended versions of the convection parametrization scheme at ECMWF and their impact on the mean and transient activity of the model in the tropics, ECMWF Tech. Memo. 206, Eur. Cent for MediumRange Weather Forecasts, Reading, UK, 1994. a, b
Ploeger, F., Konopka, P., Günther, G., Grooß, J.U., and Müller, R.: Impact of the vertical velocity scheme on modeling transport in the tropical tropopause layer, J. Geophys. Res., 115, D03301, https://doi.org/10.1029/2009JD012023, 2010. a, b
Ploeger, F., Riese, M., Haenel, F., Konopka, P., Müller, R., and Stiller, G.: Variability of stratospheric mean age of air and of the local effects of residual circulation and eddy mixing, J. Geophys. Res.Atmos., 120, 716–733, https://doi.org/10.1002/2014JD022468, 2014. a, b
Ploeger, F., Riese, M., Haenel, F., Konopka, P., Müller, R., and Stiller, G.: Variability of stratospheric mean age of air and of the local effects of residual circulation and eddy mixing, J. Geophys. Res.Atmos., 120, 716–733, https://doi.org/10.1002/2014JD022468, 2015.
Ploeger, F. and Birner, T.: Seasonal and interannual variability of lower stratospheric age of air spectra, Atmos. Chem. Phys., 16, 10195–10213, https://doi.org/10.5194/acp16101952016, 2016. a, b, c, d, e
Rayner, N. A., Parker, D. E., Horton, E. B., Folland, C. K., Alexander, L. V., Rowell, D. P., Kent, E. C., and Kaplan, A.: Global Analyses of sea surface temperatures, sea ice, and night marine air temperature since the late nineteenth century, J. Geophys.Res., 108, 4407, https://doi.org/10.1029/2002JD002670, 2003. a
Reithmeier, C. and Sausen, R.: ATTILA: Atmospheric Tracer Transport in a Lagrangian Model, Tellus B, 54, 278–299, 2002. a, b, c, d, e, f, g, h, i
Reithmeier, C., Sausen, R., and Grewe, V.: Investigating lower stratospheric model transport: Lagrangian calculations of mean age and age spectra in the GCM ECHAM4, Clim. Dynam. 30, 225–239, https://doi.org/10.1007/s0038200702941, 2008.
Riese, M., Ploeger, F., Rap, A., Vogel, B., Konopka, P., Dameris, M., and Forster, P.: Impact of uncertainties in atmospheric mixing on simulated UTLS composition and related radiative effects, J. Geophys. Res., 117, D16305, https://doi.org/10.1029/2012JD017751, 2012. a
Roeckner, E., Brokopf, R., Esch, M., Giorgetta, M., Hagemann, S., Kornblueh, L., Manzini, E., Schlese, U., and Schulzweida, U.: Sensitivity of simulated climate to horizontal and vertical resolution in the ECHAM5 atmosphere model, J. Climate, 19, 3771–3791, 2006. a, b, c, d, e, f, g
Rolinski, S., Müller, C., Heinke, J., Weindl, I., Biewald, A., Bodirsky, B. L., Bondeau, A., BoonsPrins, E. R., Bouwman, A. F., Leffelaar, P. A., te Roller, J. A., Schaphoff, S., and Thonicke, K.: Modeling vegetation and carbon dynamics of managed grasslands at the global scale with LPJmL 3.6, Geosci. Model Dev., 11, 429–451, https://doi.org/10.5194/gmd114292018, 2018.
Sander, R., Jöckel, P., Kirner, O., Kunert, A. T., Landgraf, J., and Pozzer, A.: The photolysis module JVAL14, compatible with the MESSy standard, and the JVal PreProcessor (JVPP), Geosci. Model Dev., 7, 2653–2662, https://doi.org/10.5194/gmd726532014, 2014. a
Schoeberl, M. R., Douglass, A. R., Polansky, B., Boone, C., Walker, K. A., and Bernath, P.: Estimation of stratospheric age spectrum from chemical tracers, J. Geophys. Res., 110, D21303, https://doi.org/10.1029/2005JD006125, 2005. a
Seibert, P., Krüger, B., and Frank, A.: Parametrisation of convective mixing in a Lagrangian particle dispersion model, in: Proceedings of the 5th GLOREAM Workshop, Wengen, Switzerland, 24–26 September2001, 2002. a
Stapf, B.: Lagrangesche Behandlung der Konvektion, Thesis (Diploma), FH Regensburg, Germany, 2002.
Stenke, A. and Grewe, V.: Simulation of stratospheric water vapor trends: impact on stratospheric ozone chemistry, Atmos. Chem. Phys., 5, 1257–1272, https://doi.org/10.5194/acp512572005, 2005. a
Stenke, A., Grewe, V., and Ponater, M.: Lagrangian transport of water vapor and cloud water in the ECHAM4 GCM and its impact on the cold bias, Clim. Dynam., 31, 491–506, https://doi.org/10.1007/s0038200703475. 2008. a
Stiller, G. P., von Clarmann, T., Haenel, F., Funke, B., Glatthor, N., Grabowski, U., Kellmann, S., Kiefer, M., Linden, A., Lossow, S., and LópezPuertas, M.: Observed temporal evolution of global mean age of stratospheric air for the 2002 to 2010 period, Atmos. Chem. Phys., 12, 3311–3331, https://doi.org/10.5194/acp1233112012, 2012. a, b, c
Stohl, A., Hittenberger, M., and Wotawa, G.: Validation of the Lagrangian particle dispersion model FLEXPART against largescale tracer experimant data, Atmos. Environ., 32, 4245–4264, 1998. a, b
Stohl, A., Forster, C., Frank, A., Seibert, P., and Wotawa, G.: Technical note: The Lagrangian particle dispersion model FLEXPART version 6.2, Atmos. Chem. Phys., 5, 2461–2474, https://doi.org/10.5194/acp524612005, 2005.
Tao, M.: Atmospheric Mixing in a Lagrangian Framework, Dissertation, Schriften des Forschungzentrums Jülich, Reihe Energy and Environment, Vol. 320, available at: http://elpub.bib.uniwuppertal.de/edocs/dokumente/fbc/physik/diss2016/tao/dc1612.pdf (last access: 29 April 2019), 2016. a
Tiedtke, M.: A comprehensive mass flux scheme for cumulus parameterization in largescale models, Mon. Weather Rev., 117, 1779–1800, 1989. a, b
Tost, H., Jöckel, P., and Lelieveld, J.: Influence of different convection parameterisations in a GCM, Atmos. Chem. Phys., 6, 5475–5493, https://doi.org/10.5194/acp654752006, 2006. a, b
Tost, H.: Global Modelling of Cloud, Convection and Precipitation Influences on Trace Gases and Aerosols, University of Bonn, Germany, available at: http://hss.ulb.unibonn.de/diss_online/math_nat_fak/2006/tost_holger (last access: 29 April 2019), 2006. a
Waugh, D. W. and Hall, T. M.: Age of stratospheric air: Theory, observations, and models, Rev. Geophys., 40, 1010, https://doi.org/10.1029/2000RG000101, 2002. a, b, c
Zaucker, F., Daum, P., Wetterauer, U., Beikowitz, C., Kromer, B., and Broecker, W.: Atmospheric ^{222}Rn measurements during the 1993 NARE Intensive, J. Geophys. Res., 101, 29149–29164, https://doi.org/10.1029/96JD02029, 1996. a, b
Zhang, K., Wan, H., Zhang, M., and Wang, B.: Evaluation of the atmospheric transport in a GCM using radon measurements: sensitivity to cumulus convection parameterization, Atmos. Chem. Phys., 8, 2811–2832, https://doi.org/10.5194/acp828112008, 2008. a, b, c, d, e, f
http://en.wikipedia.org/wiki/Marsaglia_polar_method(last access: 6 May 2019)
Only for uniformly distributed pseudorandom numbers, i.e. without the Marsaglia polar method.