Tropospheric mixing and parametrization of unresolved convective updrafts as implemented in the Chemical Lagrangian Model of the Stratosphere (CLaMS v2.0)

Inaccurate representation of mixing in chemistry transport models, mainly suffering from an excessive numerical diffusion, strongly influences the quantitative estimates of the stratosphere–troposphere exchange (STE). The Lagrangian view of transport offers an alternative to exploit the numerical diffusion for parametrization of the physical mixing. Here, we follow this concept and discuss how to extend the representation of tropospheric transport in the Chemical Lagrangian Model of the Stratosphere (CLaMS).

Although the current transport scheme in CLaMS (v1.0) shows a good ability to represent transport of tracers in the stably stratified stratosphere (Pommrich et al., 2014, and the references therein), there are deficiencies in the representation of the effects of convective uplift and mixing due to weak vertical stability in the troposphere. We show how the CLaMS transport scheme was modified by including additional tropospheric mixing and vertical transport due to unresolved convective updrafts by parametrizing these processes in terms of the dry and moist Brunt–Väisälä frequencies. The regions with enhanced convective updrafts in the novel CLaMS simulation covering the 2005–2008 period coincide with regions of enhanced convection as diagnosed from the satellite observations of the outgoing longwave radiation (OLR).

We analyze how well this approach improves the CLaMS representation of CO₂ in the upper troposphere and lower stratosphere, in particular the propagation of the CO₂ seasonal cycle from the planetary boundary layer (PBL) into the lower stratosphere. The CO₂ values in the PBL are specified by the CarbonTracker data set (version CT2013B), and the Comprehensive Observation Network for TRace gases by AIrLiner (CONTRAIL) observations are used to validate the model. The proposed extension of tropospheric transport increases the influence of the PBL in the middle and upper troposphere and at the same time impacts the STE. The effect on mean age away from the troposphere in the deep stratosphere is weak.

1 Introduction

Modeling of transport from a Lagrangian perspective has gained increasing popularity in the last few decades, and not only within the atmospheric community. The chance to avoid, or at least to minimize, the numerical diffusion ever present in Eulerian numerical schemes is the strongest motivation for the Lagrangian formulation of transport. Despite the obvious advantage of the Lagrangian view separating mixing from the advective part of transport, only very few Lagrangian chemical transport models (CTMs) with explicit mixing exist so far (e.g., Collins et al., 1997; Fairlie et al., 1999; Reithmeier and Sausen, 2002; McKenna et al., 2002b; Konopka et al., 2007; Wohltmann and Rex, 2009; Pugh et al., 2012).

In this paper, the 3-D version of the Chemical Lagrangian Model of the Stratosphere (CLaMS) will be used (McKenna et al., 2002b, a; Konopka et al., 2004). The novel approach of CLaMS is its parametrization of atmospheric mixing, especially in the stratosphere where vertical mixing is extremely weak due to a strong vertical stability. Whereas the common approach is to minimize the numerical diffusion ever present in the modeling of transport, CLaMS is a first attempt to apply this “undesirable disturbing effect” to parametrize the true physical mixing. This idea is realized by using some scaling properties of numerical diffusion, which are the same as those of atmospheric diffusivity (Konopka et al., 2007, 2012; Pommrich et al., 2014) and by applying numerical regridding only to the strongly deformed parts of the Lagrangian grid, where physical mixing is expected anyway (McKenna et al., 2002b; Konopka et al., 2004). This novel parametrization is an extension of the well-known pure Lagrangian, i.e., trajectory-based representation of transport (e.g., Stohl et al., 2005; Legras et al., 2005; Bowman et al., 2007; Wernli and Davies, 1997).

However, in the current version of CLaMS v1.0 (Pommrich et al., 2014) only deformations within quasi-isentropic layers driven by horizontal strain and vertical shear rates are taken into account. This is certainly a good approximation in a stably stratified stratosphere suppressing vertical mixing but not in the troposphere where vertical mixing is expected. To extend the CLaMS idea of “sufficiently strong” almost-isentropic deformations triggering mixing to the troposphere, we use the concept of atmospheric stability. The flow is called unstable if a small perturbation at initial time will exponentially grow during the course of the evolution of the flow.

One of the widely used parameters quantifying instabilities is the gradient Richardson number, Ri, describing the onset of instabilities driven by wind shear and/or buoyancy (e.g., Turner, 1973; Stull, 1988). Ri is defined as $\mathit{\text{Ri}}=N^{\mathrm{2}}/\left[\right(\mathrm{d}u/\mathrm{d}z{)}^{\mathrm{2}}+(\mathrm{d}v/\mathrm{d}z{)}^{\mathrm{2}}]$ with du∕dz and dv∕dz denoting the vertical shear of the horizontal wind components. N is the (dry) Brunt–Väisälä frequency quantifying buoyancy-driven turbulence in terms of the potential temperature lapse rate, i.e., $N^{\mathrm{2}}=(g/\mathit{\theta })\mathrm{d}\mathit{\theta }/\mathrm{d}z$ (θ – potential temperature, g – gravity of Earth, and z – geometric altitude). The flow becomes dynamically unstable or even turbulent when Ri<Ri_c (Ri_c critical value of Ri). This occurs either when the wind shear is strong enough to outweigh any stabilizing buoyant forces (denominator in the definition of Ri is large) or when the dry or, more general, moist environment is statically unstable (numerator in the definition of Ri is small or even negative because the lapse rate dθ∕dz is small or even negative). The scale-dependent value of Ri_c is about 0.25, although reported values have ranged from roughly 0.2 to 1.0 (Balsley et al., 2008).

In the stratosphere, where the flow is characterized by high static stability, only almost-isentropic deformations driven by the horizontal strain and vertical shear are considered in the current version of CLaMS v1.0 (Pommrich et al., 2014, and the references therein). These deformations measured in terms of the Lyapunov exponent λ are used in CLaMS to parametrize mixing within layers, which are parallel to the isentropes at all levels above 300 hPa. However, the effect of vertical instabilities being a dominant feature of tropospheric transport is not taken into account. To parametrize such (potentially) vertically unstable regions by using “sufficiently small” values of the dry or moist Brunt–Väisälä frequency, N or N_m, is our main heuristic idea to extend to extend the CLaMS transport scheme from the deformation-driven scheme to the deformation- and vertical-instability-driven scheme.

Generally, CLaMS tropospheric tracers like CO or CH₄ show reasonable mixing ratios in the Upper Troposphere and Lower Stratosphere (UTLS) only if their values in the orography-following lowest model boundary (thickness: 1.1–1.3 km) are artificially enhanced by a factor 1.5–2 (Pommrich et al., 2014). The consequence of this enhancement is that vertical gradients of these tracers are too strong, especially in the lower and middle troposphere if compared with observations (not shown). A possible reason might be that a significant part of the upward vertical transport within the troposphere is underestimated.

Figure 1 shows a comparison in terms of the e90 tracer between the current CLaMS version (Fig. 1a), the here discussed extension of tropospheric transport (Fig. 1b) with the Whole Atmospheric Community Climate Model (WACCM, red isolines), which is known for a good representation of tropospheric transport and chemistry (e.g., see Park et al., 2013, and the references therein). Here, the CLaMS annual means (calculated for 2007) are compared with the WACCM 1955–2099 climatology derived from the REF-C2 reference simulations, which have coupled atmosphere, ocean, ice, and land components (Abalos et al., 2017) (year-to-year variability in WACCM climatology is much smaller if compared with deviation of CLaMS from the WACCM climatology).

Figure 1 The different shades of blue quantify the annual mean of the e90 tracer mixing ratios (calculated for 2007) as a function of latitude and log-pressure altitude and are derived from the current version of CLaMS v1.0 described in Pommrich et al. (2014) (a) and from the version discussed here with extended tropospheric transport (b). For comparison, the WACCM 1955–2099 climatology of e90 is also shown (red isolines) as described in Abalos et al. (2017). The mean tropopause is depicted by the black line.

Similar to WACCM, the artificial e90 tracer, with a constant e-folding lifetime of 90 d, is set to 150 ppb everywhere in the lowest layer of CLaMS. The e90 tracer is suitable to diagnose typical timescale of transport from the planetary boundary layer (PBL) into the lower stratosphere (Prather et al., 2011; Abalos et al., 2017). Figure 1 shows that CLaMS in the current version significantly underestimates the upward transport if compared with the WACCM model and that this comparison improves if the new version of transport is included.

In this paper, we aim to parametrize the unresolved processes like convective updrafts and tropospheric mixing whose representations in global reanalysis data are uncertain (Russo et al., 2011). Although more complex convective schemes exist both in the Eulerian (e.g., Tiedtke, 1989; Emanuel, 1991) as well as in the Lagrangian formulation (e.g., Collins et al., 2002; Erukhimova and Bowman, 2006; Forster et al., 2007; Pugh et al., 2012; Ueyama et al., 2018; Brinkop and Jöckel, 2019), our approach mainly intends to cover the range of possible variability due to unresolved tropospheric transport. Guided by diabatic rather than kinematic thinking, and by the wish to reduce the complexity of convection as much as possible, we aim here to create a technical framework which allows us to estimate the impact of such uncertainties on the composition of the UTLS as well as on the Stratosphere–Troposphere exchange (STE). Even small changes in the concentration and distribution of radiatively active gases in the UTLS such as water vapor or ozone significantly impact radiative forcing on surface temperature (Riese et al., 2012).

In the next section, we describe the properties and mean distributions of the dry or moist Brunt–Väisälä frequency, N or N_m, as derived from the meteorological data. N and N_m are used in this paper to parametrize the additional tropospheric transport. Section 3 explains technical details of this parametrization. Section 4 describes the CLaMS setup (v2.0) and details the performed model simulations as well as some results diagnosing which regions of the atmosphere are mainly affected by our extension of transport. To validate these new properties, we discuss the transport of CO₂ from the PBL into the lower stratosphere and compare the respective distributions with the airborne observations and the CarbonTracker model simulations. Finally, we discuss our results in Sect. 5.

2 Vertically unstable troposphere versus stably stratified stratosphere

Static stability can be quantified in terms of the (dry) Brunt–Väisälä frequency (BVF) via $N^{\mathrm{2}}=(g/\mathit{\theta })(\mathrm{d}\mathit{\theta }/\mathrm{d}z)$. N describes the frequency at which an air parcel oscillates when displaced vertically in a statically stable environment, i.e., within a region with a positive lapse rate $\mathrm{d}\mathit{\theta }/\mathrm{d}z>\mathrm{0}$ (for some details see Appendix A). Because the well-mixed troposphere is characterized by low values of N and the stably stratified stratosphere by high values of N, it is expected that this difference also manifests in the corresponding vertical diffusivities (large and small for the troposphere and the stratosphere, respectively).

Figure 2 The zonal mean distribution of dry, N² (a, b), and moist, $N_{\mathrm{m}}^{\mathrm{2}}$ (c, d), Brunt–Väisälä frequencies (BVFs) in the troposphere and the lower stratosphere for DJF (a, c) and JJA (b, d) as calculated from 2005 ERA-Interim data. The thick yellow contours in the top row highlight those parts of the atmosphere where the tropospheric mixing is included in CLaMS. The bluish colored regions in the bottom panels mark places from which additional parametrization of convection lifts the CLaMS air parcels from the PBL into the upper troposphere.

To take into account the contribution of latent heat release to the vertical instabilities, N² can be modified by introducing the moist Brunt–Väisälä frequency, N_m, with $N_{\mathrm{m}}^{\mathrm{2}}=(g/{\mathit{\theta }}_{\mathrm{e}})(\mathrm{d}{\mathit{\theta }}_{\mathrm{e}}/\mathrm{d}z)$ where θ_e is equivalent potential temperature, i.e., the temperature an air parcel would reach if all the water vapor in the air parcel were to condense, releasing its latent heat, and then were brought down adiabatically to 1000 hPa (for details see also Appendix A). The zonal mean of the dry and moist BVF distributions during boreal winter and summer 2005 is shown as an example in Fig. 2.

Thus, lowest values of the dry BVF are distinctive for the middle troposphere especially in the tropics. The minimum in lapse rate, dθ∕dz, above the PBL represents the level of maximum convective impact (Gettelman and de Forster, 2002), which is also characterized by a vertically broad maximum of tropospheric signatures with enhanced CO and strongly reduced ozone (Pan et al., 2017, and references therein). Obviously, the values of N² clearly increase by crossing the tropopause but they also increase in the region below 700 hPa and are the lowest in the tropical and subtropical middle troposphere (i.e., within the yellow contour in Fig. 2a and b). In the following section, we use this middle tropospheric minimum of N², although zonally resolved, to parametrize tropospheric mixing in CLaMS and, in this way, to include unresolved mixing processes in the middle and upper troposphere into the CLaMS transport. Finally, we note that simple zonal means shown in Fig. 2, instead of tropopause-based zonal means, smear out the strong stability contrasts across the tropopause (Birner, 2006; Birner et al., 2006). However, we use the climatology shown in Fig. 2 only to explain our procedure; the presented transport extension is based on full 3-D stability fields which contain all the details of the tropopause as far as they are resolved by the reanalysis.

On the other hand, low values of the moist BVF (i.e., bluish regions with $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ in Fig. 2c and d) can be found at altitudes below 700 hPa between 30^∘ S and 30^∘ N, i.e., in the tropical and subtropical PBL. The air parcels with $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ are conditionally unstable (see Appendix A), i.e., they may undergo strong vertical uplift, if some other favorable condition causing saturation may happen within such air masses (e.g., gravity-wave-induced temperature fluctuations which are not resolved in the meteorological data). In the following section, we show how conditional instability can be used to trigger additional (advective) upward transport due to unresolved convection.

3 Extension of transport scheme: tropospheric mixing and unresolved convection

To extend the CLaMS mixing scheme, we follow two heuristic ideas. First, due to a much lower vertical stability in the troposphere than in the stratosphere, we enhance tropospheric mixing in the model everywhere (dry) vertical stability is sufficiently small. Second, we take into account additional transport driven by convection, especially by deep convection which is not sufficiently resolved in the reanalysis data.

Thus, whereas the first approach is related to changes in the mixing part of CLaMS and affects the nearest neighbors of each Lagrangian air parcel, the second goal is related to changes in the advection part of CLaMS, i.e., modification of the trajectory calculation. Both extensions should be driven by instabilities quantified in terms of the dry and moist BVFs, respectively, which were introduced in the previous section. By including such a revised transport scheme, we seek a better representation of transport in the free troposphere, which also likely improves the performance of the model within the UTLS region. Because all our changes are confined to the troposphere, we expect a weak influence on stratospheric transport in CLaMS which has been successfully validated in many previous studies. Furthermore, the scheme should not give a heavy burden on the computation time compared to the current version of CLaMS. Before going into the details, we shortly describe the standard version of CLaMS (in the following denoted as the reference setup).

3.1 Reference setup

As the reference, we use the 2005–2008 time slice of the 40-year CLaMS transient simulation starting on 1 January 1979 and driven by the horizontal winds and diabatic heating rates (vertical velocities) derived from the ERA-Interim reanalysis (Dee et al., 2011; Ploeger et al., 2010). This configuration as well as the model initialization follows the model setup described in Pommrich et al. (2014) (100 km horizontal and 400 m vertical resolutions around 380 K, CLaMS v1.0). The first 10 years of the CLaMS transient simulation are the model spinup. An example of the mean age of air (AoA) zonal distribution, calculated relative to the Earth's surface, is shown for 19 August 2005 as a function of the latitude and the hybrid potential temperature ζ (Fig. 7a) and will be used for comparison with CLaMS simulations including the extended tropospheric transport (Sect. 4).

The vertical coordinate ζ is the hybrid σ–θ which allows us to resolve transport processes in the troposphere influenced by the orography and transport processes in the stratosphere, where adiabatic horizontal transport dominates (Mahowald et al., 2002). More precisely, we replace η by σ in the hybrid η–θ coordinate as proposed by Mahowald et al. (2002):

$$\begin{array}{}\text{(1)}& \mathit{\zeta }:=f\left(\mathit{\sigma }\right)\mathit{\theta }(p,T),\end{array}$$

with

$$\begin{array}{}\text{(2)}& f\left(\mathit{\sigma }\right)=\left\{\begin{array}{rlrlr}{\displaystyle }& {\displaystyle }\sin \left({\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}-\mathit{\sigma }}{\mathrm{1}-{\mathit{\sigma }}_{\mathrm{r}}}}\right)& {\displaystyle }\mathit{\sigma }& {\displaystyle }>{\mathit{\sigma }}_{\mathrm{r}},& {\displaystyle }\\ {\displaystyle }& {\displaystyle }\mathrm{1}& {\displaystyle }\mathit{\sigma }& {\displaystyle }\le {\mathit{\sigma }}_{\mathrm{r}},& {\displaystyle }{\mathit{\sigma }}_{\mathrm{r}}={\displaystyle \frac{p_{\mathrm{r}}}{p_{\mathrm{0}}}}.\end{array}\right.\end{array}$$

Here, p₀ denotes the constant reference pressure level set to 1000 hPa. p_r defines the pressure level around which the dry potential temperature θ smoothly transforms into the terrain-following coordinate $\mathit{\sigma }=p/p_{\mathrm{s}}$.

For p_r, Mahowald et al. (2002) used the value 300 hPa, i.e., σ_r=0.3. For situations with no orography (e.g., see surface with p_s=1000 hPa), this means that the condition $\mathit{\sigma }=p/p_{\mathrm{s}}<{\mathit{\sigma }}_{\mathrm{r}}$ is valid everywhere above the level of 300 hPa. Consequently, in this region the vertical coordinate ζ is given by the dry potential temperature θ. Conversely, below the level of 300 hPa, θ smoothly transforms into σ. For situations with orography (e.g., at the summit of Mount Everest with p_s≈330 hPa), the condition $\mathit{\sigma }=p/p_{\mathrm{s}}<{\mathit{\sigma }}_{\mathrm{r}}$ is first valid everywhere above ≈100 hPa. Note that all p_r values between 0 and 1000 hPa are possible with the consequence that higher p_r values extend the applicability of θ as a vertical coordinate down to the lower troposphere. In the setup used here, ζ covers the range between ζ=0 (i.e., orography-following) and $\mathit{\zeta }=\mathit{\theta }=\mathrm{2500}$ K (Pommrich et al., 2014).

Thus, because ζ in the troposphere is a less intuitive coordinate (hybrid mixture between σ and θ), isolines of pressure and potential temperature as well as the zonal mean of the World Meteorological Organization (WMO) tropopause are also shown in Fig. 7a. Furthermore, vertical boundaries of the layers Δζ_i, $i=\mathrm{1},\mathrm{\dots },N$, within which CLaMS mixing is organized (for more details see Konopka et al., 2007), are depicted in the right side of Fig. 7a. In the CLaMS reference run, the lowest layer Δζ₀ approximating the PBL extends between ζ=0 (Earth's surface) and 100 K, i.e., $\mathrm{\Delta }{\mathit{\zeta }}_{\mathrm{0}}=\mathrm{\Delta }{\mathit{\zeta }}_{\mathrm{pbl}}=\mathrm{100}$ K. After each trajectory step Δt, mixing ratios of all air parcels within this layer are replaced by their initial configuration and prescribed by a lower boundary condition (the same procedure is applied for the upper boundary, Δζ_N=300 K, covering the ζ-range between 2200 and 2500 K.

In the default mixing scheme, CLaMS uses the integral deformations γ=λΔt derived from the relative motion of the nearest neighbors within each layer Δζ_i (adaptive grid procedure). Here, λ is the Lyapunov exponent of such a deformation and Δt denotes the advective time step (typically between 6 and 24 h). In the stratosphere and in large parts of the UTLS, where the flow is characterized by a high static stability, only sufficiently strong deformations with γ>γ_c (γ_c denoting an empirical critical deformation) are expected to trigger mixing with the best choice for γ_c between 0.8 and 1.5 (Konopka et al., 2004). This also means that there is some freedom in the choice of the parameters λ and Δt. Whereas for the stratosphere the values of Δt=24 h and λ=1.5 d⁻¹ (γ_c=1.5) were used in the past (Konopka et al., 2004) (and are used here as the reference configuration), a larger mixing frequency with Δt=6 h seems to give better results in the well-mixed troposphere and in the UTLS region (Vogel et al., 2011; Konopka and Pan, 2012).

3.2 Tropospheric mixing

In the following, we assume that the additional tropospheric mixing should be triggered whenever the corresponding value of N² interpolated at the CLaMS air parcel is less than a critical value denoted in the following as $N_{\mathrm{c}}^{\mathrm{2}}$. $N_{\mathrm{c}}^{\mathrm{2}}$ is a free parameter which, basically, can be adjusted by comparison with the experimental data. We expect that $N_{\mathrm{c}}^{\mathrm{2}}$ should be around zero and should identify regions with (statistically) enhanced tropospheric mixing. If an air parcel fulfills the criterion $N^{\mathrm{2}}<N_{\mathrm{c}}^{\mathrm{2}}$, the air parcel will be mixed with all nearest neighbors diagnosed by Delaunay triangulation in the respective CLaMS layer under consideration (so we use the same nearest neighbors as in the standard CLaMS mixing scheme). In this way the composition of the considered air parcels are affected, but not their geometric positions.

Figure 3 Mixing driven by vertical instability and strong wind shear. The profile on the left side is an idealized N² profile derived from the reanalysis and interpolated on the CLaMS air parcels (red). The lower part is representative for the troposphere with $N^{\mathrm{2}}<N_{\mathrm{c}}^{\mathrm{2}}$ and the upper part for the stratosphere with $N^{\mathrm{2}}>N_{\mathrm{c}}^{\mathrm{2}}$. Thus, in the lower layer static stability is weak with the opposite configuration in the upper layer where, in addition, the wind shear is strong. In the upper layer, adaptive regridding is used to include mixing between a subset of the nearest neighbors of the considered air parcel (default mixing in CLaMS). In addition, in the lower layer, all nearest neighbors will be mixed with the considered air parcel if criterion $N^{\mathrm{2}}<N_{\mathrm{c}}^{\mathrm{2}}$ is valid. The purple parcels are mixed parcels of red and blue parcels.

Figure 3 is a schematic diagram illustrating how the tropospheric extension of mixing works. If this additional mixing is applied, we set the new mixing ratio of a considered air parcel and its nearest neighbors through averaging their composition, which is shown as a change of the parcel's color. This setting completes the mixing without changing any model parcel position and the number of parcels because it can be executed directly in the current mixing module (more precisely, after the deformation-driven adaptive grid procedure of CLaMS).

3.3 Unresolved (deep) convection

Commonly, convection is understood as a vigorous vertical updraft resulting from instabilities in the vertical temperature and water distribution profiles, preferably connecting the PBL with the free troposphere and is schematically shown in Fig. 4a.

Figure 4 (a) Schematic definition of convective fluxes and scales. Vertically unstable profile of the equivalent potential temperature θ_e (black thin line) triggers convective updraft in our simplified scheme. (b) Convective updraft in CLaMS: the original trajectory (thick solid red) and the modified trajectory (dashed red) with convective transport from the lowest CLaMS layer Δζ_pbl at time step t₀+Δt. The vertical displacement Δθ of the considered air parcel is estimated through the latent heat release of condensation and may happen at every time step along the trajectory when the reanalysis data are read (i.e., every Δt=6 h for ERA-Interim, see text for more details).

Convection starts on relatively small horizontal scales on the order of a few kilometers and is mainly driven by the latent heat release of the gaseous and liquid water content. Such diabatic processes may lift a large amount of air. The respective upward mass flux is traditionally denoted as an updraft (thick upward red arrow in Fig. 4a).

The reverse diabatic mass transport, the so-called downdraft, is related to the evaporation of water or melting of ice and is in general much smaller than the original convective updraft (thin blue arrows). A much larger part of the downward mass flux, the so-called subsidence, occurs on much larger horizontal scales of the order of 10–100 km, or even larger if Rossby waves or gravity waves are induced by convection (thick downward red arrows). A significant part of subsidence results mainly from longwave cooling of radiatively active water vapor which is vertically transported and freeze-dried within the convective tower and subsequently spread horizontally at the level of the main convective outflow. Note that the horizontal scale of this outflow region is at least 1 order of magnitude larger than the horizontal extension of the convective tower at the ground.

It is generally believed that the exchange of mass driven by deep convection can efficiently inject the air masses from the PBL into the upper troposphere or even, although very rarely, into the lower stratosphere (Schiller et al., 2009; Corti et al., 2006). In fact, the extension of mixing presented in the previous subsection is still limited by the model layers and, consequently, is not suitable to parametrize unresolved convective events which connect the PBL with layers in the upper troposphere on a timescale of minutes to hours.

We now present an alternative method to enhance upward transport for conditionally unstable air parcels with $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ in order to lift such air masses from the lowest layer of CLaMS Δζ_pbl (following the orography and approximating here the PBL) into the upper troposphere. Figure 4b shows the concept of estimating the uplift of boundary air by adding Δθ to the trajectory in the vertical direction when the condition $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ is diagnosed along the trajectory (i.e., using diabatic concept of vertical motion).

Following Ertel (1938) (for details see Appendix B), we use the following approximation for Δθ:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }\mathit{\theta }={\displaystyle \frac{L_{\mathrm{v}}{\mathit{\theta }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{w}}}{c_{p}T}},\end{array}$$

where θ₀ and μ_w denote the potential temperature and the total water-vapor mass mixing ratio in the air parcel where the condition $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ is fulfilled. T denotes the temperature in the air parcel and here approximates the saturation temperature. L_v is the specific latent heat for evaporation and c_p denotes the specific heat at constant pressure, p.

Figure 5 Color coded is the fraction of ERA-Interim time steps (6 h frequency) within a season when the criterion triggering the deep convection scheme is fulfilled at CLaMS air parcels within the lowest layer of the model Δζ_pbl=250 K (σ_r=0.7, for details see text). Left and right columns show DJF and JJA 2005 climatologies, respectively. In the top row the isolines of the outgoing longwave radiation (OLR) as derived from the NOAA satellite archive are highlighted while in the bottom row the isolines of the convective available potential energy (CAPE) are overlaid (cyan).

To illustrate how such a parametrization works, the zonally resolved fraction of events with $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ occurring within the lowest CLaMS layer are shown in Fig. 5. The respective DJF and JJA climatologies derived from ERA-Interim for 2005 reveal the expected spatial distribution, land-ocean contrasts, and the seasonality. To justify the use of $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ as a proxy of convection, we compare its spatial distribution with the satellite-based outgoing longwave radiation (OLR, Fig. 5a and b, cyan isolines) as well as with the ERA-Interim-based convective available potential energy (CAPE) as an alternative method to detect convection (for definition of CAPE, see Appendix B). The comparison shows a good correlation between the climatology of the regions with $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ and the respective OLR (Fig. 5a and b), whereas the correlation with the CAPE is less pronounced (e.g., in the region around 30–50^∘ E, 45^∘ N during JJA). Motivated by this finding and because of the simplicity of tracing conditionally unstable air parcels, we use in the following the criterion $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ as the first condition to trigger convective events.

Our second condition is related to the question of whether every conditionally unstable air parcel is a source of convection which should be taken into account. This question is also related to the fact that the number of air parcels in CLaMS is not strictly conserved but kept roughly constant within about ±10 % flexibility through the adaptive regridding procedure (current mixing scheme). It means that the mixing procedure is able to adjust a certain increase or decrease in the number of air parcels, but this amount should be below ±10 %. Figure 6 shows the 1-year climatological PDF of Δθ in the tropics (30^∘ S–30^∘ N). The probability of Δθ larger than 35 K is around 30 % and it decreases rapidly from 35 to 60 K. When the Δθ is too small to leave the lower boundary, it is not necessary to add Δθ to the trajectory. Thus, as our second condition, we only uplift such an air parcel along the trajectory (i.e., within the advective step) if Δθ is sufficiently large. In our control runs, the minimum of Δθ is set to Δθ_crit=35 K and is shown as the dashed line in Fig. 6. This choice is also the reason that we call this parametrization a “deep” convection scheme.

Figure 6 The probability distribution function (PDF, in units of % K⁻¹) of convective uplift Δθ in the tropics (30^∘ S-30^∘ N) for the air parcels in the lowest level of CLaMS (Δζ=250 K, Δz≈1.4–2 km) where $N_{\mathrm{m}}^{\mathrm{2}}$ is negative using the ERA-Interim reanalysis data of the year 2005. The restriction for the deep convection scheme ($\mathrm{\Delta }\mathit{\theta }>\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{crit}}=\mathrm{35}$ K) is marked as the dashed straight line.

To estimate the influence of our convection parametrization on the mass budget, we calculate the annually and globally averaged mass fluxes due to parametrized convective updrafts (for details see Appendix C). Because only updrafts are parametrized, this has an impact on the total mass budget in CLaMS with a surplus of the annually and globally averaged mass flux density of around 10 kg m⁻² d⁻¹ between the upper edge of the lowest CLaMS layer and around 360 K (highest level of the here parametrized convection). On the other hand, the used ERA-Interim diabatic vertical velocities are per se not mass conserving. Fueglistaler et al. (2009) have shown that there is a significant deficit of mass on the order of 20 kg m⁻² d⁻¹ (cf. their Fig. 10), roughly in the same altitude range where the CLaMS convective parametrization works.

A possible explanation for this deficit can be related to the fact that the horizontal resolution of the ERA-Interim reanalysis (roughly 80 km) does not sufficiently resolve the convective towers (which are of the order of 1 km) but does better resolve the large-scale subsidence which may be better reproduced in the ERA-Interim diabatic vertical velocities. Thus, our parametrization aims to close this imbalance by including some additional convective updrafts. The results show that qualitatively our simple approach roughly balances such a deficit and there is still some potential to enhance the strength of parametrized convection. At least the influence of our parametrization on the mass budged is comparable with “intrinsic uncertainties” in the mass budget of the reanalysis itself. Following the procedure described in Rosenlof (1995), and applied in Konopka et al. (2010) (see their appendix), it is possible to restore the mass balance violations at least in the annual mean. It is also worth noting that unresolved convection is a process not present in the divergence of the horizontal wind and, consequently, not influencing the 3-D (large-scale) velocity field.

Finally, it should be mentioned that our Lagrangian parcels are still too large (around 100 km in horizontal and a few hundred meters in vertical direction) to be transported by realistic convective systems. Thus, they are more suitable to describe large-scale convective outflow rather than the convective towers which are well below their horizontal resolution. This fact justifies to a certain extent our restriction to consider only deep convective events with Δθ>35 K. Note also that by increasing model resolution by a factor of 2, both horizontally and vertically, the respective convective mass flux would decrease by a factor of 8.

4 CLaMS performance with additional tropospheric transport

In this section, we describe the details of the CLaMS configuration for the simulations with extended tropospheric transport (v2.0), show in which part of the atmosphere the CLaMS air parcels are affected by this extension, and compare the respective AoA distributions. We especially investigate the propagation of the CO₂ distribution from the boundary layer into the lower stratosphere for different model configurations and evaluate the related CO₂ variability (annual cycle and trend) with the observations and the CarbonTracker model simulations.

4.1 Setup for reference and control simulations

As for the reference simulation described in subsection 3.1 (in the following denoted as REF), all other simulations cover the same time period 2005–2008 and have the same vertical and horizontal resolutions above the tropopause. All these simulations start on 1 January 2005 and are initialized with the REF concentrations. The duration of 3 years for all our simulation is long enough to see the changes in the troposphere and tropical lower stratosphere. However, 3 years is certainly too short for the whole stratosphere, but this is our compromise concerning the available computing time and is long enough to see at least in which direction the changes of the model have to be expected. Table 1 provides the key information for all CLaMS simulations discussed in this paper.

Figure 7 Zonal mean of mean age of air (AoA) for the reference simulation (a) and for the simulations with the full extension of tropospheric transport (b, c) calculated for 19 August 2005 (i.e., after more than 8 months of transport). For the definition of the abbreviations see Table 1. The hybrid vertical coordinate ζ is used and the isolines of the potential temperature θ (black), pressure p (cyan) and the zonal mean of the WMO tropopause are also shown. On the right side of each panel, the boundaries of the CLaMS layers, Δζ_i, $i=\mathrm{1},\mathrm{\dots },N$, are plotted (see text for more details).

Table 1 List of CLaMS reference and control simulations with different configurations of mixing, unresolved convective updrafts, and parameter σ_r defining the transition between the σ- and θ-related vertical coordinate. For all simulations with unresolved convective updrafts and tropospheric mixing, the critical values of $N_{\mathrm{m}}^{\mathrm{2}}$ and N² are set to 0 and $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻², respectively, and Δθ_crit=35 K.

In addition to the REF configuration, we use two slightly different configurations: σ_r=0.7, Δζ_pbl=250 K (REF_0.7, Fig. 7b) and σ_r=0.3, Δζ_pbl=140 K (REF_0.3, Fig. 7c) with the letter choice being very close to the REF case with σ_r=0.3 and Δζ_pbl=100 K (see below). Grid configurations for all cases listed in Table 1 are the same in the stratosphere although significant differences are in the troposphere: by using σ_r=0.7, isentropic mixing and diabatic vertical velocities, two central concepts of CLaMS, can be extended to a larger part of the troposphere, in particular to the middle tropical troposphere. Note that the 320 and 330 K isentropes in Fig. 7b are within CLaMS layers defined by the ζ coordinate and, consequently, mixing within such layers is roughly isentropic. Thus, almost the whole UTLS region, down to the tropical middle troposphere, is covered by such isentropic layers in the σ_r=0.7 simulation (see Tao et al., 2018, especially their Appendix 1).

Furthermore, both apparently different choices of Δζ_pbl for REF_0.7 and REF_0.3 correspond roughly to the same geometric thickness of the lowest CLaMS layer, which varies between 1.4 and 2.2 km and approximates here the PBL. Because in the REF case this thickness is certainly too small (between 1.1 and 1.3 km), we decided to enhance these values to more realistic numbers. Thus, for a fairer comparison, we take REF_0.7 or REF_0.3 as the main reference (instead of REF) with σ_r=0.7 being our main choice. It should be emphasized that by using predefined boundary conditions in the PBL, we do not resolve any transport in this part of the atmosphere and confine our efforts only to improve transport in the free troposphere extending between the PBL and the tropopause.

As a first step toward extended tropospheric transport (in the following denoted as control runs), we decrease the advective time step of trajectories from 24 to 6 h and, to keep the intensity of the standard CLaMS mixing scheme roughly constant, we also increase the Lyapunov exponent from 1.5 to 3.5 d⁻¹ corresponding to the critical deformation γ_c=0.85 (case MIX_0.7 in Table 1). Such a slightly higher mixing frequency relative to the reference case (every 6 h instead of every 24 h) has proven to give better representation of the CO–ozone correlations in the UTLS region (Vogel et al., 2011; Konopka and Pan, 2012). Higher mixing frequency also resolves the diurnal cycle of our new parametrization of tropospheric transport described in the last section.

The next step is to add new tropospheric transport, i.e., tropospheric mixing (TROP_MIX) and unresolved convective updrafts (UNRES_CONV). Before discussing these new contributions we first show our results in Fig. 7.

Here, the zonal mean distribution of AoA, calculated relative to the Earth's surface for 1 day (19 August 2005), is plotted for the three cases discussed above: REF (Fig. 7a) and two control simulations with full extension of tropospheric transport, FULL_EXT_0.7 (Fig. 7b) and FULL_EXT_0.3 (Fig. 7c). Both types of control simulations show much younger air in the troposphere if compared with the reference run. Also the gradients across the tropopause are more pronounced. On the other hand, stratospheric distributions are very similar for all three cases. Now we go more into the details of such a tropospheric extension.

4.2 Diagnostic of extended tropospheric transport

It is easy to tag and count all air parcels in CLaMS which undergo additional tropospheric mixing and which are lifted from the lower boundary to the middle and upper troposphere by the deep-convection scheme introduced in the previous section. In Fig. 8a and b the DJF and JJA zonally averaged fractions of additionally mixed air parcels (calculated for 2005) are color coded as a function of latitude and pressure. The black isolines in Fig. 8a and b approximate the fraction of the CLaMS air parcels which was lifted from the lowest boundary layer to the middle and upper tropical troposphere using the deep convection parametrization. For comparison, the fractions of the number of air parcels affected by the default CLaMS mixing scheme are shown in Fig. 8c and d, which always happens, independent of whether the discussed extension of transport is included or not. In addition, the respective mean WMO tropopause as well as the isentropes is shown.

Figure 8 (a, b) DJF/JJA 1-year (2005) climatology of percentage of CLaMS air parcels undergoing tropospheric mixing (colors) and of air parcels lifted from the lowest layer of the model into the middle and upper (tropical) troposphere (black contours). (c, d) Same type of climatology but for air parcels which undergo the standard CLaMS mixing procedure (adaptive regridding driven by horizontal strain and vertical shear within the CLaMS layers). The respective WMO tropopause (beige), horizontal wind marking the position of the jets (light gray), and the isentropes (dark gray) are also shown.

Note that CLaMS tropospheric mixing practically does not affect any air parcels in the stratosphere (the zero line of the calculated fraction is well below the tropopause, not shown). Note also that numbers of air parcels affected by the deep convection scheme are smaller than 10 % with highest levels (around 360 K) during JJA, mainly related to the Asian summer monsoon (not shown). Furthermore, both tropospheric mixing and the deep convection transport show some seasonality like N² and $N_{\mathrm{m}}^{\mathrm{2}}$, respectively (i.e., higher in the summer hemisphere). Finally, the default mixing scheme is much weaker than the tropospheric mixing although the seasonality is very similar. Note that this part of mixing is also present in the vicinity of the tropopause although with stronger signatures on the tropical side of the jets and below the tropopause. At least in this climatological picture, only standard mixing in CLaMS contributes to a direct STE.

4.3 Validation with CO₂ observations

CO₂ is a useful tracer for validation of transport in the models, mainly in the troposphere and lower stratosphere where CO₂ is basically chemically inert (Waugh and Hall, 2002). The only stratospheric source of CO₂ is a small contribution (<1 ppmv) from methane oxidation (Bates and Nicolet, 1950; Ehhalt, 1974; Boucher et al., 2009) that is taken into account in all CLaMS simulations discussed here (Pommrich et al., 2014). Thus, the quality of the CO₂ distribution reproduced in CLaMS is determined largely by the quality of the lower boundary condition and the quality of the representation of transport. For the latter, tropospheric transport is a significant part.

The atmospheric mixing ratios of CO₂ are essentially both monotonically increasing (trend) and periodic (seasonality) signals which define a stringent test for the model representation of tropospheric transport and STE (Bönisch et al., 2008, 2009). As recently shown by Diallo et al. (2017), even inverted vertical profiles of CO₂ across the extratropical tropopause are possible during the northern hemispheric summer despite a continuous increase in the mean CO₂ in the PBL resulting from the growing anthropogenic emissions.

In CLaMS, CO₂ mixing ratios propagate upwards from the lowest layer, Δζ_pbl, for which the CarbonTracker data set was used (Peters et al., 2007; the updates are documented at http://carbontracker.noaa.gov, last access: 20 June 2019) with CO₂ mixing ratios available for the 2000–2012 period (simulation CT2013B, available every 3 h, see ftp://products/carbontracker/co2/CT2013B/molefractions/co2_total/, last access: 20 June 2019). In particular, the first five lowest levels of each CarbonTracker data set were vertically averaged and used to overwrite CLaMS air parcels within the PBL layer every 6 h. The reference simulation was initialized at 1 January 2000 and beginning from 1 January 2005 all other control simulations were started using the output of the reference simulation for the initial distribution.

The zonal means of CO₂, also calculated for two representative days, 5 May and 25 September, 2005, are shown in Fig. 9. In particular, results for the reference simulation (REF) and for the two control simulations with full tropospheric transport (FULL_EXT, ${\mathit{\sigma }}_{\mathrm{r}}=\mathrm{0.3}/\mathrm{0.7}$) can be compared with the respective CarbonTracker distribution which was used in CLaMS to initialize the lower boundary of the model. In all CO₂ distributions, the upward propagation of the annual cycle can be clearly diagnosed with higher values during the boreal summer and vertical inversion during the fall. However, the propagation of the tropospheric signal shows some obvious differences with a faster upward propagation in CLaMS control simulations than in the CLaMS reference run. Note that the cross-hemispheric transport is weaker in the σ_r=0.7 than in the σ_r=0.3 CLaMS configuration. Note also that the upward propagation of the CO₂ annual cycle is well-confined by the position of the tropopause (black dots) in all CLaMS runs, while in the CarbonTracker data this property is less pronounced although, as in CLaMS, ERA-Interim reanalysis is used in the underlying transport model (see CT2013B documentation: https://www.esrl.noaa.gov/gmd/ccgg/carbontracker/CT2013B, last access: 20 June 2019). We will come to this point later.

Figure 9 Upward propagation of the CO₂ distribution from the lowest layer of the model where CO₂ was initialized by the CarbonTracker data (CT2013B) on two days: 1 May 2005 (a–d) and 25 September 2005 (e–h). The CLaMS zonal means for different model configurations are compared with the CarbonTracker distribution itself. Black dots denote the tropopause derived from the ERA-Interim data.

Now, the CO₂ time and space evolution derived from CLaMS simulations are compared with the observations of the Comprehensive Observation Network for TRace gases by AIrLiner (CONTRAIL) (Machida et al., 2008). CO₂ mixing ratios were measured during regular flights by Japan Airlines from Japan to Australia, Europe, North America, and Asia with continuous measuring equipment (CME) for in situ CO₂ observations, as well as improved automatic air sampling equipment (ASE) for flask sampling (for more details about the instruments see Machida et al., 2002). This data set provides significant spatial coverage, particularly in the Northern Hemisphere (Sawa et al., 2015). CONTRAIL observations have a vertical resolution of a few meters (during ascents and descents) and a horizontal resolution of a few hundred meters, resulting from the high sampling frequency of these instruments.

Figure 10 (a, b) Seasonality of CO₂ for CLaMS simulations compared with CONTRAIL observations (black filled circles and vertical lines denote mean values and their standard deviations, respectively) at two locations: 15^∘N and 5.5 km (a) and 15^∘N and 10.5 km (b). (c) The accumulated errors in CO₂ for CLaMS simulations compared with all CONTRAIL observations as a function of the altitude.

Here, we use the zonally and monthly averaged time evolution of these observations between 2005 and 2008 interpolated at a latitude–altitude grid with 10^∘ × 1 km resolutions and extending between 20^∘ S and 60^∘ N and between 5.5 and 12.5 km, respectively (for more details see Diallo et al., 2017). Comparison of these mean CONTRAIL observations with the respective CLaMS results for the reference and all control simulations are shown in Fig. 10.

In particular, the comparison of the seasonal cycle and trend at 15^∘ N for two selected altitudes of 5.5 and 10.5 km is plotted in Fig. 10a and b. The right panel shows the accumulated error, i.e., the zonal mean of the mean square deviation between the CLaMS simulation and CONTRAIL observations averaged over all latitude–altitude grid points where gridded (mean) measurements and their standard deviations are available. Note a higher variability in CO₂ at 5.5 km compared to 10.5 km, indicating a stronger variability due to the differences in the sources of CO₂ than due to the vicinity of the tropopause where our gridding procedure does not differ between the stratospheric and tropospheric values. Thus, the reference simulation (REF) and the simulation with 6 h mixing frequency (MIX_0.7) show not only an amplitude that is too small but also that their phase is delayed if compared with the CONTRAIL observations.

There is a clear improvement in the representation of the CO₂ distribution quantified in terms of the phase and the amplitude of the seasonal cycle as well as in terms of the accumulated error by taking additional tropospheric transport into account. The best results are achieved by including both the tropospheric mixing and the convection parametrization (orange and red curves are for FULL_EXT_0.7 and FULL_EXT_0.3, respectively). By switching off the tropospheric mixing or the convection parametrization, or both (here the results only for σ_r=0.7 are shown), the cumulative error increases up to 80 %. While the tropospheric mixing is more important below 9 km, the improvement due to convection parametrization dominates between 9 and 13 km. Although there are still some errors in the amplitude, the additional tropospheric transport significantly improves the overall agreement. Note that also the CarbonTracker distribution, even the one achieved by assimilating observations, does not show a perfect comparison with the CONTRAIL observations (which are not included in the assimilation procedure of the version used here, CT2013B). Remarkably, CLaMS control simulations are becoming even better than CarbonTracker distributions in the region above 10 km (cf. Fig. 9) probably caused by a very limited vertical resolution of the CarbonTracker data around the tropopause (only 6 levels between 9 and 18 km).

4.4 Impacts on the stratosphere

For this purpose, we discuss in Fig. 11 the differences in the distribution of AoA due to extension of tropospheric transport by considering its annual and zonal mean calculated for the year 2007 (last year of our simulations covering the 2005–2007 period).

Figure 11 Annual and zonal mean of AoA (for 2007, in years) as derived from the reference simulation REF_0.7 (a) and the differences relative to this distributions (in month) calculated for some configurations listed in Table 1. Black line denotes the tropopause.

As a reference distribution we use the AoA of the REF_0.7 case (shown in the Fig. 11a) and calculate the respective differences in the AoA distributions of the simulations with 6 h default mixing (MIX_0.7, Fig. 11b), with added unresolved convective updrafts (UNRES_CONV_0.7, Fig. 11c), with added tropospheric mixing (TROP_MIX_0.7, Fig. 11d), with the full tropospheric extension (FULL_EXT_0.7, Fig. 11e) and with full tropospheric extension but slightly reduced default mixing by increasing λ_c from 3.5 to 4.0 d⁻¹ (FULL_EXT_SM_0.7, Fig. 11f).

As expected, the air below the tropical tropopause becomes younger by up to 6 months if the additional tropospheric transport is included (Fig. 11c to f). The effect of convection is mainly confined to the tropics while the strongest impact of tropospheric mixing can be diagnosed in the high latitudes. However, and at first surprisingly, the air becomes slightly older in the stratosphere (around 6 months). Note that this is not a consequence of enhanced tropospheric transport but of the change in the default mixing scheme from Δt=24 h and λ_c=1.5 d⁻¹ (reference) to Δt=6 h and λ_c=3.5 d⁻¹ and can be deduced from Fig. 11b. Such a change leads to a slightly enhanced isentropic mixing across the tropical pipe, which also enhances the stratospheric recirculation and makes the stratospheric air older due to aging by mixing (Garny et al., 2014; Poshyvailo et al., 2018). Consistently, aging by mixing in the stratosphere becomes smaller by reducing the isentropic part of CLaMS mixing by setting Δt=6 h and λ_c=4.0 d⁻¹ in the default mixing scheme (Fig. 11f) although the tropospheric AoA is only weakly affected. This indicates a secondary role of CLaMS standard mixing scheme in the troposphere and underlines its primary role in the stratosphere. In this way, the cross-tropopause gradient of AoA can be enhanced by combining the stratospheric mixing with the proposed extension of tropospheric transport.

5 Conclusions

Implementation of mixing in Lagrangian transport models is still an important issue in the ongoing scientific discussion. Here, we follow the idea of using numerical diffusion to parametrize physical mixing which was first proposed and implemented in connection with the Chemical Lagrangian Model of the Stratosphere (CLaMS). In particular, we extend this idea to the troposphere where vertical stability is much smaller if compared with the stratosphere for which CLaMS was originally developed. By using the lapse rates of the dry and moist potential temperature mainly defining the squares of the dry and moist Brunt–Väisälä frequencies N² and $N_{\mathrm{m}}^{\mathrm{2}}$, we parametrize two important tropospheric processes which are not sufficiently resolved in the current version of CLaMS v1.0 (Pommrich et al., 2014), i.e., tropospheric mixing in regions with small lapse rates of the dry potential temperature and unresolved (deep) convection in regions with conditionally unstable lapse rates of the moist potential temperature.

The implementation of both processes improves CLaMS performance measured here in terms of the quality of the simulated CO₂ seasonality. However, there is still some freedom in the choice of the free parameters (like critical values for N², $N_{\mathrm{m}}^{\mathrm{2}}$, or Δθ) and in the representation of the PBL in the model. In particular the choice of the critical displacement Δθ triggering the onset of convective updrafts (currently set to 35K) needs further investigations. Using smaller values for this parameter would enhance the violation of the mass conservation and would make it mandatory to include some diabatic correction to the vertical winds following the procedure described in Rosenlof (1995) and Konopka et al. (2010).

By including other species like CO, ozone, CH₄, and water vapor, and comparing such distributions with observations, we plan in the future to reduce these degrees of freedom. On the other hand, by covering the whole range of possible variability of tropospheric transport, we will also be able to find the respective variability in the air composition entering the stratosphere.

Appendix A: Vertical instability, Brunt–Väisälä frequency, and convective available potential energy (CAPE)

Vertical instability is strongly related to the concept of buoyancy. We now condense some textbook knowledge and start from the Boussinesq approximation of the vertical momentum equation by taking into account only buoyancy effects (e.g., see Salby, 1996, or Vallis, 2006):

$$\begin{array}{}\text{(A1)}& {\displaystyle \frac{Dw}{Dt}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\partial \mathrm{\Delta }p}{\partial z}}-g{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}}\approx -g{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}}.\end{array}$$

The solution of this equation describes the vertical velocity w of an air parcel in a hydrostatic reference atmosphere, i.e., defined by the relation ${\partial }_z{p}_{\mathrm{0}}=-{\mathit{\rho }}_{\mathrm{0}}g$ with mean pressure and air mass density profiles given by p₀(z) and ρ₀(z), respectively (z – geometric altitude, g – gravity of Earth). The Δ quantities, i.e., Δp and Δρ, describe the deviation from such a mean state.

The right-hand side of Eq. (A1) quantifies the buoyancy. Following Archimedes, the mass-related buoyancy F_b∕m can be derived from the weight of the environmental fluid replaced by the parcel:

$$\begin{array}{}\text{(A2)}& \begin{array}{rl}-{\displaystyle \frac{F_{\mathrm{b}}}{m}}& ={\displaystyle \frac{\mathit{\rho }gV-{\mathit{\rho }}_{\mathrm{0}}gV}{\mathit{\rho }V}}=g{\displaystyle \frac{\mathit{\rho }-{\mathit{\rho }}_{\mathrm{0}}}{\mathit{\rho }}}\stackrel{\mathrm{Boussinesq}}{=}g{\displaystyle \frac{\mathit{\rho }-{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}}\\ & =g{\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}},\end{array}\end{array}$$

where, using the Boussinesq approximation, ρ was replaced by ρ₀. With $w=\mathrm{\Delta }z/\mathrm{d}t$ and $\mathrm{\Delta }\mathit{\rho }=-(\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}/\mathrm{d}z)\mathrm{\Delta }z$, the following equation for Δz can be derived:

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{{\mathrm{d}}^{\mathrm{2}}}{\mathrm{d}{t}^{\mathrm{2}}}}\mathrm{\Delta }z+N^{\mathrm{2}}\mathrm{\Delta }z=\mathrm{0},N^{\mathrm{2}}=-{\displaystyle \frac{g}{{\mathit{\rho }}_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{d}z}},\end{array}$$

with N denoting the Brunt–Väisälä frequency. Using the ideal gas law and assuming the incompressibility of the flow, we get (for details see Vallis, 2006)

$$\begin{array}{}\text{(A4)}& {\displaystyle \frac{\mathrm{\Delta }\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{0}}}}=-{\displaystyle \frac{\mathrm{\Delta }\mathit{\theta }}{{\mathit{\theta }}_{\mathrm{0}}}}.\end{array}$$

With $\mathrm{\Delta }\mathit{\rho }=-\mathrm{\Delta }\mathit{\theta }({\mathit{\rho }}_{\mathrm{0}}/{\mathit{\theta }}_{\mathrm{0}})=-({\mathit{\rho }}_{\mathrm{0}}/{\mathit{\theta }}_{\mathrm{0}})(\mathrm{d}{\mathit{\theta }}_{\mathrm{0}}/\mathrm{d}z)\mathrm{\Delta }z$, the Brunt–Väisälä frequency can also be rewritten with the well-known definition (index “0” is omitted):

$$\begin{array}{}\text{(A5)}& N^{\mathrm{2}}={\displaystyle \frac{g}{\mathit{\theta }}}{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}z}}.\end{array}$$

In general, Eq. (A3) has two types of solutions: periodic (i.e., stable) solutions for N²>0 and exponentially increasing (i.e., unstable) solutions for N²<0. We conclude that the lapse rate of the potential temperature, θ, defines the stability and the instability of the atmospheric environment through the positive and negative lapse rate dθ∕dz (or through the positive and negative values of N²), respectively.

Figure A1 The equivalent potential temperature T_e of an air parcel with temperature T and pressure p and with the dry mass m_d, water vapor content m_v, and the liquid and/or solid water content m_s. T_e is defined by a reversible process completely removing m_v by vapor-water condensation (i.e., transforming pseudoadiabatically m_v into m_vs) and by using this energy to heat the original air parcel by p=const, i.e., T_e>T. Liquid and/or solid water (m_s) are removed, but without any correction of T_e. (i.e., due to sedimentation this is a pseudoadiabatic and not a pure adiabatic process).

Now we generalize this concept to the atmosphere containing water vapor, i.e., to the moist atmosphere (e.g., see Salby, 1996). First, we define the equivalent potential temperature θ_e by using the equivalent temperature T_e, i.e., the temperature of an air parcel from which all the water vapor has been extracted by an adiabatic vapor-water condensation process (see also Fig. A1):

$$\begin{array}{}\text{(A6)}& {\mathit{\theta }}_{\mathrm{e}}={\displaystyle \frac{T_{\mathrm{e}}}{T}}\mathit{\theta },T_{\mathrm{e}}=T+{\displaystyle \frac{L_{\mathrm{v}}}{c_p}}{\mathit{\mu }}_{\mathrm{w}}.\end{array}$$

L_v is the latent heat of evaporation and μ_w the water vapor mixing ratio. The energy released by the phase transition from liquid water to the ice phase can be neglected in most cases because the respective latent heat is smaller by a factor of 10 than the latent heat of gas–liquid transition (334 kJ for melting versus 2270 kJ for evaporation for 1 kg liquid water). There is a number of different definitions of the equivalent potential temperature (Bolton, 1980). Our definition corresponds to the simplified formula proposed by Stull (1988).

Using the same type of arguments as for the dry atmosphere, we also quantify the vertical instability of the moist atmosphere in terms of the lapse rate of the equivalent potential temperature θ_e or in terms of the respective (moist) Brunt–Väisälä frequency N_m

$$\begin{array}{}\text{(A7)}& N_{\mathrm{m}}^{\mathrm{2}}={\displaystyle \frac{g}{{\mathit{\theta }}_{\mathrm{e}}}}{\displaystyle \frac{\mathrm{d}{\mathit{\theta }}_{\mathrm{e}}}{\mathrm{d}z}}.\end{array}$$

However, atmospheric environments with a negative lapse rate of the equivalent potential temperature or with negative values of $N_{\mathrm{m}}^{\mathrm{2}}$ define only the so-called conditionally unstable atmosphere (e.g., see Salby, 1996), i.e., regions which could be unstable if the respective phase transition releasing latent heat would happen (such air parcels with $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ are not necessarily saturated, so some unresolved motions like adiabatic gravity waves are needed to get saturation).

For comparison, we also use the known concept of convective available potential energy (CAPE), which can be understood as a different measure of the unstable buoyancy (Emanuel, 1994). Starting form Eqs. (A2) and (A4), we can write

$$\begin{array}{}\text{(A8)}& -{\displaystyle \frac{F_{\mathrm{b}}}{m}}=g{\displaystyle \frac{\mathit{\rho }-{\mathit{\rho }}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}}}=-g{\displaystyle \frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{0}}}}.\end{array}$$

Figure A2 CAPE (red) versus stability-based (cyan) measure of potential convective uplift, Δθ, transporting CLaMS air parcels from the PBL to the upper troposphere. For more explanation see the text.

Then, CAPE is defined as the following integral (in J kg⁻¹, see also Fig. A2):

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}\mathrm{CAPE}& =\underset{z_{min}}{\overset{z_{max}}{\int }}\text{“unstable buoyancy force” }\mathrm{d}z\\ & =\underset{\mathrm{LFC}}{\overset{\mathrm{EL}}{\int }}g{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{wb}}\left(z\right)-{\mathit{\theta }}_{\mathrm{0}}\left(z\right)}{{\mathit{\theta }}_{\mathrm{0}}\left(z\right)}}\mathrm{d}z.\end{array}\end{array}$$

Whereas θ₀(z) is the dry potential temperature of the environment (ambient air), θ_wb(z) is the so-called wet-bulb potential temperature which needs some further explanations: we begin with the lifting condensation level (LCL) defined as the height at which a parcel of air becomes saturated when it is lifted adiabatically from the Earth's surface (so the potential temperature does not change). Starting from the LCL the air parcel is then transported along a moist adiabat (also known as saturation-adiabatic process, i.e., a pseudoadiabatic process for which the air is saturated). The corresponding dry potential temperature along such a moist adiabat defines the wet-bulb potential temperature. Note that a distinction is made between the reversible process, in which total water is conserved, and the pseudoadiabatic or irreversible moist adiabatic process, in which liquid water is assumed to be removed as soon as it is condensed (see also Fig. A1). The cross points of such a moist adiabate with θ₀(z) define the region with unstable conditions (red). The lowest cross point defines the level of free convection (LFC), whereas the equilibrium level (EL) is the height where the (potential) temperature of a buoyantly rising parcel again equals the (potential) temperature of the environment.

In this paper, we use the condition $N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$ rather than condition CAPE>0 to find those air parcels in the lowest model level (approximating the PBL) which might undergo convection (see Fig. A2). In the next appendix, we discuss an additional assumption quantifying the vertical displacement of such air parcels.

Appendix B: Latent heat release versus vertical transport

To derive Eq. (3), quantifying the convection-driven vertical uplift Δθ from the latent heat δQ available within the air parcel, we assume that the source of heating is only the latent heat release of water vapor condensation. For a unit mass, we get along a moist adiabate (e.g., see Salby, 1996)

$$\begin{array}{}\text{(B1)}& \mathit{\delta }Q=-L_{\mathrm{v}}\mathrm{d}{\mathit{\mu }}_{\mathrm{s}},\end{array}$$

where L_v is the specific latent heat for water evaporation (or condensation, in J kg⁻¹) and μ_s is the saturation water-vapor mass mixing ratio.

Using the entropy definition of the potential temperature $\mathit{\theta }:={\mathit{\theta }}_{\mathrm{0}}\exp \left(s/c_p\right)$ with s being the specific entropy measured in J (K kg)⁻¹ (i.e., entropy per unit mass), c_p denoting the specific heat at constant pressure, p, and θ₀ being the temperature, T₀, that the air parcel would acquire if adiabatically brought to the surface pressure, p₀, we get

$$\begin{array}{}\text{(B2)}& \mathrm{d}s=c_p{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathit{\theta }}}.\end{array}$$

Thus, using the second law of thermodynamics, $\mathrm{d}s=\mathit{\delta }Q/T$, and substituting δQ by the Eq. (B1), dθ can be derived as

$$\begin{array}{}\text{(B3)}& \mathrm{d}\mathit{\theta }=-{\displaystyle \frac{L_{\mathrm{v}}\mathit{\theta }\mathrm{d}{\mathit{\mu }}_{\mathrm{s}}}{c_{p}T}},\end{array}$$

where T is the saturation temperature. Furthermore, we assume that during the condensation process, when the saturated air parcel moves from the initial state θ₀ to the final state θ₀+Δθ, the saturation temperature T does not change. Thus, Eq. (B3) can be simply integrated

$$\begin{array}{}\text{(B4)}& \underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{0}}+\mathrm{\Delta }\mathit{\theta }}{\int }}{\displaystyle \frac{\mathrm{d}\mathit{\theta }}{\mathit{\theta }}}=-\underset{{\mathit{\mu }}_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{0}}\right)}{\overset{{\mathit{\mu }}_{\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}}+\mathrm{\Delta }\mathit{\theta })}{\int }}{\displaystyle \frac{L_{\mathrm{v}}\mathrm{d}{\mathit{\mu }}_{\mathrm{s}}}{c_{p}T}}.\end{array}$$

The result of the integral can be written as

$$\begin{array}{}\text{(B5)}& \mathrm{\Delta }\mathit{\theta }={\mathit{\theta }}_{\mathrm{0}}\left[\exp \left({\displaystyle \frac{L_{\mathrm{v}}\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{s}}}{c_{p}T}}\right)-\mathrm{1}\right],\end{array}$$

where $\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{s}}={\mathit{\mu }}_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{0}}\right)-{\mathit{\mu }}_{\mathrm{s}}({\mathit{\theta }}_{\mathrm{0}}+\mathrm{\Delta }\mathit{\theta })$. Because $\frac{L_{\mathrm{v}}\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{s}}}{c_{p}T}\ll \mathrm{1}$, a first-order approximation of the right-hand side of Eq. (B5) through the Taylor expansion of the exponential function (Ertel, 1938) is

$$\begin{array}{}\text{(B6)}& \mathrm{\Delta }\mathit{\theta }={\displaystyle \frac{L_{\mathrm{v}}{\mathit{\theta }}_{\mathrm{0}}\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{s}}}{c_{p}T}},\end{array}$$

which gives a relationship between the total change in the potential temperature and the change of water-vapor mass mixing ratio. Strictly, Δμ_s is the change of the water vapor saturation mass mixing ratio before and after a model time step. In our deep convection scheme, Δμ_s is estimated by the total water-vapor mass mixing ratio, μ_w, before the model time step with the assumption that (1) the timescale of deep convection and its associated condensation is smaller than one model time step (here 6 h) and (2) the residual water vapor content after deep convection (μ_s(θ₀+Δθ)) is so small that it can be neglected. Thus, Δμ_s is assumed to be the total water-vapor mass mixing ratio μ_w in the air parcel within the lowest layer of CLaMS where the criterion ($N_{\mathrm{m}}^{\mathrm{2}}<\mathrm{0}$) is fulfilled. Therefore, the uplifting of an air parcel in our deep convection scheme is estimated by

$$\begin{array}{}\text{(B7)}& \mathrm{\Delta }\mathit{\theta }={\displaystyle \frac{L_{\mathrm{v}}{\mathit{\theta }}_{\mathrm{0}}{\mathit{\mu }}_{\mathrm{w}}}{c_{p}T}}.\end{array}$$

Appendix C: Mass flux due to unresolved convective updrafts

To estimate the influence of our convection parametrization on the mass budget, we calculate the zonally resolved (i.e., per latitude bin Δϕ) mass density change, ${\dot{M}}_{\mathrm{c}}=\mathrm{d}{M}_{\mathrm{c}}/\mathrm{d}t$, within CLaMS layer Δζ for all convective events which were initiated during one advection time step $\mathrm{\Delta }t=\mathrm{1}/\mathrm{4}$ d:

$$\begin{array}{}\text{(C1)}& {\dot{M}}_{\mathrm{c}}(\mathit{\varphi },\mathit{\zeta })=\mathit{\sigma }(\mathit{\varphi },\mathit{\zeta })f_m(\mathit{\varphi },\mathit{\zeta }){\displaystyle \frac{\mathrm{\Delta }\mathit{\zeta }}{\mathrm{\Delta }t}}.\end{array}$$

Units of ${\dot{M}}_{\mathrm{c}}$ are kg m⁻² d⁻¹. The ζ-related air density σ is given as $-(\mathrm{1}/g)\partial p/\partial \mathit{\zeta }$ (units: kg km⁻²). Furthermore, f_m is the fraction of air parcels convectively transported within the latitude bin Δϕ to a target layer Δζ denoted here as n_conv relative to the total number of air parcels in this destination region denoted here as N_tot, i.e., $f_m=n_{\mathrm{conv}}/N_{\mathrm{tot}}$ (we assume that all air parcels are uniformly distributed within the layer, both horizontally and vertically). The related convective updraft mass flux density, F_c, is given through ${\partial }_z{F}_{\mathrm{c}}={\dot{M}}_{\mathrm{c}}$ or as a discretized solution:

$$\begin{array}{}\text{(C2)}& F_{\mathrm{c}}({\mathit{\varphi }}_i,{\mathit{\zeta }}_j)=\sum _{k=j}^{k=N}{\dot{M}}_{\mathrm{c}}({\mathit{\varphi }}_i,{\mathit{\zeta }}_k),\end{array}$$

where N denotes the total number of CLaMS layers and ϕ_i, $i=\mathrm{1},\mathrm{\dots },M$, are the latitude bins.

Figure C1 shows the annually averaged values of F_c (Fig. C1a) and its global mean (Fig. C1b). The latter profile quantifies the additional mass flux into the model caused by parametrized convective updrafts and it roughly balances the mass deficit diagnosed in Fueglistaler et al. (2009) (cf. their Fig. 10) due to use of the non-mass-conserving diabatic vertical velocities.

Figure C1 Annually averaged and zonally resolved updrafts due to parametrized convection (σ_r=0.3) (a) and their global mean (b). On the right side of panel (a) the boundaries of CLaMS layers, Δζ_i, i=1…N, are shown.