ATTILA 4.0: Lagrangian advective and convective transport of passive tracers within the ECHAM5/MESSy (2.53.0) chemistry–climate model

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 grid-point submodel LGTMIX to calculate the mixing of compounds in Lagrangian representation. The new infrastructure routines simplify the transformation between grid-point (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 ²²²Rn on Lagrangian parcels. To evaluate the new physical routines, two simulations in free-running 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 age-of-air distribution is more in accordance with observations. The global tropospheric distribution of ²²²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 inter-parcel mixing on the age of air in the tropopause region and the stratosphere, but there is no significant effect for the troposphere.

1 Introduction

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 large-scale 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 3-D 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 climate-optimised 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 large-scale 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 grid-scale temperature and humidity and provides a matrix for vertical convective particle displacement (Seibert et al., 2002; Forster et al., 2007).

In the former (non-parallelised) version of ATTILA, convective tracer tendencies were calculated in grid-point 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 ²²²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 Model description

2.1 EMAC – a MESSy-fied 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 multi-institutional computer codes. The core atmospheric model is the fifth-generation 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 T42L47MA-resolution, 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 set-up 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)

Table 1 List of MESSy submodels used for the simulations in this study.

^* Formerly called OFFLEM.

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 ²²²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 sub-submodel that hosts the basic transformation routines needed in ATTILA to convert variables from grid-point to Lagrangian representation and vice versa (ATTILA_TOOLS),

  • –

    a new submodel that uses ATTILA_TOOLS to calculate the transformation of user-specified variables between Lagrangian and grid-point 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 flux-form semi-Lagrangian (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 (time-stepping 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 message-passing 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):

$$\begin{array}{}\text{(1)}& N=\sum _{i=\mathrm{0}}^{p-\mathrm{1}}{n}_i,\end{array}$$

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 Lagrangian-enabled submodels as so-called 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 MPI-based 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 grid-point representation. The latter is solely used to transform the Lagrangian tracers into grid-point space for output and further analyses.

Subroutines to transform and transpose variables between Lagrangian representation and (parallel decomposed) grid-point representation, and vice versa, are collected in a specific toolbox module named ATTILA_TOOLS. This also comprises specific subroutines for the transformation of grid-point 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.

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

2. 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).

3. 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).

4. 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 pseudo-random 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 pseudo-random number series is used for the initial distribution of the parcels in the model atmosphere. These pseudo-random number series are provided by the MESSy infrastructure submodel RND. RND provides uniformly distributed pseudo-random 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¹. The generation of high-quality pseudo-random 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 pseudo-random numbers to all other tasks. This implies a load imbalance and requires additional MPI communication, yet for most pseudo-random number generators it is the only possibility. However, for the Mersenne Twister² (among others) Haramoto et al. (2008) found an efficient “jump ahead” facility, i.e. a method to advance the pseudo-random 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 pseudo-random numbers. Jumping ahead by numbers not representable in the form a×2^b can be achieved by additionally harvesting r pseudo-random 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 pseudo-random number series and has been implemented in the MESSy infrastructure submodel RND. The same pseudo-random 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 pseudo-random numbers it needs to harvest. The only prerequisite for this to work is that the number of required pseudo-random numbers per (each) task is a priori known to all other tasks. For instance, if for each ATTILA parcel k pseudo-random numbers are required (e.g. per model time step), in total k×N pseudo-random numbers need to be harvested, i.e. k×n_i for task i. That means that task i needs to jump ahead by

$$\begin{array}{}\text{(2)}& j_i^{\prime }=k\times \sum _{q=\mathrm{0}}^{i-\mathrm{1}}{n}_i\end{array}$$

steps before it can harvest its own chunk of k×n_i pseudo-random numbers.

For the simulations analysed below, we used three uniformly distributed pseudo-random 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 three-dimensional wind field using a fourth-order 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,

$$\begin{array}{}\text{(3)}& \mathit{\xi }=\mathit{\theta }f,\end{array}$$

with θ being the potential temperature and f being defined as

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{If}p>p_{\mathrm{r}}f=\sin \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}-\frac{p}{p_{\mathrm{s}}}}{\mathrm{1}-\frac{p_{\mathrm{r}}}{p_{\mathrm{0}}}}}\right),\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{If}p\le p_{\mathrm{r}}f=\mathrm{1},\end{array}$$

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₀=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):

$$\begin{array}{}\text{(6)}& \dot{\mathit{\xi }}=\dot{\mathit{\theta }}f+\mathit{\theta }\dot{f}.\end{array}$$

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 grid-box 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 mid-level 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^uX^u, F^d=M^dX^d, and F^e=M^eX^e:

$$\begin{array}{}\text{(7)}& -{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}{{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{\mathit{\rho }}\stackrel{\mathrm{‾}}{w^{\prime }{X}^{\prime }}\right)}{\partial z}}}_{\mathrm{conv}}=-{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{‾}}{\mathit{\rho }}}}\left({\displaystyle \frac{\partial F^{\mathrm{u}}}{\partial z}}+{\displaystyle \frac{\partial F^{\mathrm{d}}}{\partial z}}+{\displaystyle \frac{\partial F^{\mathrm{e}}}{\partial z}}\right),\end{array}$$

where $\stackrel{\mathrm{‾}}{X}$ is the tracer mass mixing ratio, M is the mass flux, $\stackrel{\mathrm{‾}}{\mathit{\rho }}$ is the air density, $\stackrel{\mathrm{‾}}{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 grid-box 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.

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle \frac{\partial M^{\mathrm{u}}}{\partial z}}=E^{\mathrm{u}}-D^{\mathrm{u}}\text{(9)}& {\displaystyle }& {\displaystyle \frac{\partial M^{\mathrm{d}}}{\partial z}}=E^{\mathrm{d}}-D^{\mathrm{d}}\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{with}{M}^{\mathrm{e}}=-(M^{\mathrm{u}}+M^{\mathrm{d}})\end{array}$$

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 large-scale moisture convergence in cases of deep or mid-level 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.

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle \frac{\partial F^{\mathrm{u}}}{\partial z}}={\displaystyle \frac{\partial M^{\mathrm{u}}{X}^{\mathrm{u}}}{\partial z}}=E^{\mathrm{u}}{X}^{\mathrm{e}}-D^{\mathrm{u}}{X}^{\mathrm{u}}\text{(12)}& {\displaystyle }& {\displaystyle \frac{\partial F^{\mathrm{d}}}{\partial z}}={\displaystyle \frac{\partial M^{\mathrm{d}}{X}^{\mathrm{d}}}{\partial z}}=E^{\mathrm{d}}{X}^{\mathrm{e}}-D^{\mathrm{d}}{X}^{\mathrm{d}}\text{(13)}& {\displaystyle }& {\displaystyle \frac{\partial F^{\mathrm{e}}}{\partial z}}={\displaystyle \frac{\partial \left(M^{\mathrm{e}}{X}^{\mathrm{e}}\right)}{\partial z}}\end{array}$$

The calculation of the convective transport of tracer mass starts with the determination of the type of convection (deep, shallow, mid-level). 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

$$\begin{array}{}\text{(14)}& P_k={\displaystyle \frac{m_e}{m_g}}={\displaystyle \frac{(M_k-M_{k+\mathrm{1}})\mathrm{\Delta }tg}{p_{k+\mathrm{1}}-p_k}},\end{array}$$

with

$$\begin{array}{}\text{(15)}& m_g={\displaystyle \frac{(p_{k+\mathrm{1}}-p_k)A}{g}}\text{and}{m}_e=(M_k-M_{k+\mathrm{1}})A\mathrm{\Delta }t,\end{array}$$

where p is pressure (hPa), A is area (m²), M represents air mass fluxes ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\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:

$$\begin{array}{}\text{(16)}& P=(M_k-M_{k+\mathrm{1}})/M_{k+\mathrm{1}}.\end{array}$$

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 large-scale 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.

Figure 1 Mode of operation of Lagrangian convection in a vertical column. Coloured circles are Lagrangian parcels; n=3, 4, and 5 is the original number of parcels in a grid box (chosen arbitrarily for this example) that should be reached again after the convective event to keep the air mass in each grid box constant.

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 grid-point 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 grid-point 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 parcel-to-parcel communication, we define a background mixing ratio $\stackrel{\mathrm{‾}}{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:

$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{‾}}{c}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{c}_{\mathrm{i}}.\end{array}$$

The altered mixing ratio of the respective parcel is then calculated by $c_{\mathrm{i}}^{\mathrm{new}}=c_{\mathrm{i}}+(\stackrel{\mathrm{‾}}{c}-c_{\mathrm{i}})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 2-D iso-surfaces in 3-D scalar fields and to map 3-D scalar fields in GP representation on iso-surfaces; 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 grid-point 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 two-dimensional surface (e.g. isentropic surface, potential vorticity iso-surface, 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:

$$\begin{array}{}\text{(18)}& F_{\mathrm{sfc}}\left(c\right)={\displaystyle \frac{\sum m_{\mathrm{i}}\times c_{\mathrm{i}}}{\mathrm{\Delta }tA}},\end{array}$$

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 ²²²Rn and its decay products as tracers in Lagrangian representation (see Sect. 2.2.1). Here, we simulate ²²²Rn with a constant ²²²Rn source of 10 000 $\mathrm{atoms}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ over ice-free land (zero elsewhere) and a decay with a half-life of 3.8 d as the only ²²²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 ²²²Rn mass into parcels which reside lowest in the boundary layer. An example &CPL namelist is shown in the Supplement.

3 Observations

3.1 ²²²Rn

²²²Rn has been frequently used for an evaluation of large-scale 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 ²²²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). ²²²Rn is emitted from land surfaces due to a radioactive decay of radium in soils. ²²²Rn has a characteristically short radioactive half-life (τ_½=3.8 d). For the evaluation of the simulated ²²²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 ²²²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 ²²²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₆) (Stiller et al., 2012; Haenel et al., 2015) and CO₂ (Engel et al., 2009; Andrews et al., 2001). Both will be used for comparison.

4 Model simulations

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 RC1-base-08 free-running 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₂, O₃, N₂O, CF₂Cl₂, CFCl₃) have been prescribed from RC1-base-08. Methane was initialised as in RC1-base-08 and the same pseudo-emission 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¹D have been prescribed as monthly averages from RC1-base-08. The photolysis rate J(CH₄) was calculated with the submodel JVAL.

ATTILA was initialised with 1.15×10⁶ 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 grid-point 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.

• ²²²Rn is a commonly used tracer to study vertical transport into the upper troposphere due to its characteristically short radioactive half-life of (τ_½=3.8 d). In our simulations it is emitted at the surface with an emission rate of 10 000 $\mathrm{atoms}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ over ice-free land surfaces between 90^∘ N and 90^∘ S.

• SF₆_AoA and SF₆_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₆_AoA the linear-in-time increasing mixing ratio is latitude dependent. Using a spatially constant surface source (SF₆_AoAc) has the advantage that concentration differences in the atmosphere cannot have their origin in the distribution of surface sources.

• SF₆_AoA_nm has the same properties as the tracer SF₆_AoA; however, the inter-parcel mixing was set to zero (see Sect. 2.4). Hence, it can be used in comparison to SF₆_AoA to study the influence of local mixing between adjacent parcels on the global AoA distribution.

5 Evaluation

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 ²²²Rn and SF₆_AoA tracer distributions with observations and with EMAC results, i.e. from the GP space.

5.1 Simulation of ²²²Rn

Jöckel et al. (2010) showed that the EMAC model (version 2.40) is able to realistically simulate the ²²²Rn distribution. We therefore assume here that our GP simulation with the EMAC model (now version 2.53) simulates ²²²Rn similarly. A comparison with observations will follow in the next section. The inter-comparison of our new simulations of ²²²Rn between GP and LG space (Fig. 2) shows the largest values of ²²²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 ²²²Rn emissions. In the boundary layer (north of 40^∘ N) the zonal mean ²²²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 ²²²Rn values are relatively small. At 100 hPa in Fig. 2 the largest ²²²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 ²²²Rn has largely decayed due to its short half-life and the relatively long transport times in the stratosphere.

Figure 2 Zonal mean ²²²Rn mixing ratio in the boundary layer (averaged over the lowest three model layers; dashed line, left vertical axis) and at 100 hPa (solid line, right vertical axis). GP denotes the grid-point simulation with EMAC (black line), LG(diab) the simulation with EMAC–ATTILA using the diabatic vertical velocity (red line), and LG(kin) the EMAC–ATTILA simulation with kinematic vertical velocity (light blue line).

Figure 3 Zonal mean ²²²Rn mixing ratio between 800 and 1013 hPa (weighted by the level thickness) over land (dashed) and over sea (solid) for GP (black line), LG(diab) (red line), and LG(kin) (light blue line).

5.1.1 Annual cycle of ²²²Rn at the surface layer

We use ²²²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 long-term averages of ²²²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 long-range transport (²²²Rn lower than 1000 mBq m⁻³; 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 ²²²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 set-up with 90 vertical levels, whereas in our study the simulation is free running and our model set-up 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 ²²²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 ²²²Rn values for Kerguelen are not captured in the simulation. The simulated ²²²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.

Figure 4 Monthly mean ²²²Rn mBq m⁻³(STP) surface concentrations: GP and LG model results against observations at different sites from Zhang et al. (2008). The thick dashed lines include a range within a factor of 2 of the observations.

5.1.2 Vertical profiles of ²²²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 ²²²Rn in the atmosphere for several flights. The thick dashed curve is a spline interpolation on the grid of all flights. The ²²²Rn mixing ratio decreases with height up to about 3 km and remains around 10⁻²⁴ mol mol⁻¹. 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 ²²²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 ²²²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.

Figure 5 Vertical profiles of ²²²Rn mixing ratio (10⁻²¹ mol mol⁻¹) during the NARE campaign. The thick dashed line (limited to a height of 6000 m) is a spline interpolation of the scattered measurement data (triangles) on the grid. Thin dashed lines represent the one σ standard deviation of the respective simulations.

Figure 6 Vertical profiles of ²²²Rn mixing ratio (10⁻²¹ mol mol⁻¹) during the Moffett campaign. The thick dashed line is a spline interpolation of the scattered measurement data (circles) on the grid.

5.2 Age of air

The calculation of AoA is performed in two different ways: by a so-called clock tracer (a linear-in-time increasing tracer like SF₆) 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 inter-parcel 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 clock-tracer distribution can be directly compared to AoA from an observed tracer distribution.

Figure 7 Zonal mean AoA at 20 km of height (∼ 50 hPa) from the SF₆_AoA tracer (red: LG(diab), green: GP, blue: LG(kin)) as a mean over the years 2000–2010. Thick black line with circles: MIPAS data as a mean over the years 2003, 2007, 2008, 2009, 2010, and 2011. Triangles: AoA from SF₆ at 20 km; from Waugh and Hall (2002). Thin black dashed lines: minimum and maximum AoA from CO₂ at 20 km; from Andrews et al. (2001).

5.2.1 Mean age of air

Mean AoA in the stratosphere is calculated from the SF₆_AoA tracer. The transit time is estimated by comparing the tracer mixing ratio in the stratosphere with the SF₆_AoA mixing ratio in the tropical boundary layer. Figure 7 shows a comparison of SF₆_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₆ tracer, as well as from CO₂ 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₆ 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 (1960-2010) of LG(diab) is shown in Fig. 8. It confirms the well-known 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 flux-form semi-Lagrangian 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).

Figure 8 Zonal mean AoA (in years) from the SF₆_AoA tracer in the LG(diab) simulation.

Figure 9 Zonal mean difference of AoA (LG(diab)–GP) from the SF₆_AoAc tracer (with a spatially constant surface source of SF₆). The non-stippled area is statistically significant to the 99 % level.

Figure 10 Normalised age spectra (tropics and poles) of the years 2006–2010 between 50 and 0.1 hPa from the LG(diab) simulation. The data points represent the normalised frequency per half-year bin.

Figure 11 Normalised seasonal age spectrum from LG(diab) simulation of the years 1990–2010 between 400 and 500 K at 50–70^∘ N.

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 so-called 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.5-year 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.

Figure 12 Zonal mean difference between LG(diab) with standard mixing of the SF₆_AoA tracer and LG(diab) with no mixing (nm) between adjacent parcels. The non-stippled area is statistically significant to the 99 % level.

Figure 13 Zonal mean difference between LG(diab) and LG(diab-nm) with no inter-parcel mixing. The stippled area is statistically significant to the 99 % level.

5.2.3 Sensitivity of age of air to inter-parcel mixing

AoA is influenced by the amount of mixing between adjacent parcels. Inter-parcel mixing can be regarded as a diffusion process leading to a reduction of local AoA gradients. The effect of inter-parcel mixing makes stratospheric air generally younger (Fig. 12). Stratospheric AoA without inter-parcel 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 inter-parcel 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). Inter-parcel mixing in the troposphere (Fig. 13) has only a small and statistically insignificant effect on the simulated AoA because the troposphere is a well-mixed 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 inter-parcel mixing can be changed.

Figure 14 Zonal mean mass fluxes through the 380 K (black), 400 K (red), and 500 K (green) isentropes. Upward (solid line) and downward (dashed line) mass flux from the LG(diab) simulation as a mean over 1960–2010.

Figure 15 Monthly net mass fluxes from the LG(diab) simulation through the 380 K isentrope for the northern and southern extratropics.

5.3 Stratosphere–troposphere exchange

Stratosphere–troposphere exchange (STE) is characterised by a global-scale 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⁹ kg s⁻¹ and falls in the range given by Appenzeller et al. (1996) of 6–7 ×10⁹ kg s⁻¹.

Figure 16 Movement characteristic of parcels transported in the updraft in 1997. The vertical axis describes the start level of a parcel, and the horizontal axis describes the respective final updraft (end) level. Displayed are the respective numbers of parcels transported in the updraft, normalised with the maximum number.

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 set-up.

Figure 17 Similar to Fig. 16 but for parcels transported in the downdraft in 1997.

Figure 18 Similar to Fig. 16 but for parcels representing the compensating movement in the environment (subsidence) in 1997. Note that there is also a compensating upward parcel movement (shown below the blue diagonal line). For details, see the text.

6 Summary and outlook

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 inter-parcel 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 ²²²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 ²²²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 large-scale water vapour transport consistently in convective regions and to assess the convective contribution to the stratospheric water vapour budget.