How robust are stratospheric age of air trends from different reanalyses?

An accelerating Brewer–Dobson circulation (BDC) is a robust signal of climate change in model predictions but has been questioned by trace gas observations. We analyse the stratospheric mean age of air and the full age spectrum as measures for the BDC and its trend. Age of air is calculated using the Chemical Lagrangian Model of the Stratosphere (CLaMS) driven by ERA-Interim, JRA-55 and MERRA-2 reanalysis data to assess the robustness of the representation of the BDC in current generation meteorological reanalyses. We find that the climatological mean age significantly depends on the reanalysis, with JRA-55 showing the youngest and MERRA-2 the oldest mean age. Consideration of the age spectrum indicates that the older air for MERRA-2 is related to a stronger spectrum tail, which is likely associated with weaker tropical upwelling and stronger recirculation. Seasonality of stratospheric transport is robustly represented in reanalyses, with similar mean age variations and age spectrum peaks. Long-term changes from 1989 to 2015 turn out to be similar for the reanalyses with mainly decreasing mean age accompanied by a shift of the age spectrum peak towards shorter transit times, resembling the forced response in climate model simulations to increasing greenhouse gas concentrations. For the shorter periods, 1989–2001 and 2002–2015, the age of air changes are less robust. Only ERA-Interim shows the hemispheric dipole pattern in age changes from 2002 to 2015 as viewed by recent satellite observations. Consequently, the representation of decadal variability of the BDC in current generation reanalyses appears less robust and is a major uncertainty of modelling the BDC.

1 Introduction

The global circulation of the stratosphere, known as the Brewer–Dobson circulation (Brewer, 1949; Holton et al., 1995), is a crucial factor controlling the lower stratospheric composition of radiatively active trace gases; therefore, it plays an important role in the Earth's radiation budget and in climate. Hence, a realistic representation of the BDC is a prerequisite for reliable climate model predictions. However, current climate models and observations disagree regarding long-term changes of the BDC (e.g. Waugh, 2009). Climate models simulate a strengthening and accelerating circulation which is not evident from observations, representing a major uncertainty in current model predictions (see Butchart, 2014, for a recent review). The BDC is characterised by the slow upwelling in the tropics from the troposphere across the tropical tropopause layer (TTL) into the stratosphere, followed by poleward motion in the stratosphere and downwelling at middle and high latitudes. From the perspective of tracer transport, the BDC includes a residual mean mass circulation and additional two-way eddy mixing, which causes net transport of tracers but not mass (e.g. Plumb, 2002). The residual mean mass circulation may be separated into two different main branches (Birner and Bönisch, 2011). On the one hand, a shallow circulation branch transports air masses to middle latitude regions at lower levels above the tropopause within a few months to about 2 years. A deep circulation branch, on the other hand, causes transport deep into the stratosphere and downwelling at high latitudes on a timescale of several years (Birner and Bönisch, 2011). The stratospheric BDC is mechanically driven by the transfer of momentum from breaking atmospheric waves to the zonal flow (Haynes et al., 1991; Holton et al., 1995). The shallow circulation branch is mainly driven by synoptic and planetary-scale waves, and the deep branch is mainly driven by planetary waves propagating deep into the stratosphere (e.g. Plumb, 2002). Results from an idealised model suggest that the strength of the shallow circulation branch is largely controlled by the strength of wave sources in the troposphere (Gerber, 2012). The strength of the deep circulation branch, in comparison, appears to be largely controlled by the propagation conditions for waves in the stratosphere.

As the BDC is a zonal mean circulation and related wind velocities are small residuals, directly measuring the circulation is not possible (e.g. Butchart, 2014). A common measure of the BDC is the age of air, which represents the timescale for transport through the stratosphere. By definition, age of air measures the speed of the circulation. Strictly, due to atmospheric mixing processes on a multitude of scales, a stratospheric air parcel is not characterised by a single transport timescale but by a transit time distribution. This transit time distribution is termed the “age of air spectrum” (Hall and Plumb, 1994). Commonly used reference surfaces for measuring the transit time are the tropical tropopause or the tropical surface (e.g. Waugh and Hall, 2002). The first moment of the age spectrum distribution, termed the “mean age of air”, can be estimated from observations of specific trace gas species with linearly increasing mixing ratios in the troposphere, as is approximately the case for CO₂ and SF₆ (Hall and Plumb, 1994). Similar to a tracer, the mean age is affected by both the residual mean mass circulation and atmospheric mixing processes (e.g. Neu and Plumb, 1999; Garny et al., 2014), which complicates an unambiguous interpretation of mean age changes in terms of processes. The full age spectrum allows for a clearer interpretation, but its deduction from observations is much more challenging and usually hinges on simplifying assumptions like the stationarity of the flow or a specific parameterization of the age spectrum shape (e.g. Andrews et al., 1999; Schoeberl et al., 2005). Ongoing modelling studies suggest promising new methods for deducing the age spectrum based on either an improved parametric approach (Hauck et al., 2019) or an inversion approach (Podglajen and Ploeger, 2019). However, applicability to measurement data remains to be shown.

The BDC is characterised by variability on very different timescales. Seasonal variability in the BDC is caused by seasonality in wave driving, which is stronger in boreal than austral winter due to larger wave excitation by orography in the Northern Hemisphere (NH). As a result, tropical upwelling and related extratropical downwelling (in the NH) maximise in boreal winter (e.g. Yulaeva et al., 1994). On inter-annual timescales, the BDC is significantly modulated by the stratospheric Quasi-Biennial Oscillation (QBO) (e.g. Baldwin et al., 2001), by the El Niño–Southern Oscillation (ENSO) (e.g. Calvo et al., 2010; Diallo et al., 2018) and by stratospheric aerosol injected by volcanic eruptions (e.g. Diallo et al., 2017). In the long term, climate models simulate a robust strengthening of the BDC with global warming, resulting in an increase of the tropical upwelling mass flux as well as in a global decrease of the stratospheric mean age of air over the last decades and into the future (e.g. Butchart et al., 2010; Butchart, 2014). In contrast, estimates of the mean age of air from balloon observations of long-lived trace gas species in the NH middle latitudes above about 24 km show no significant long-term trend (Engel et al., 2009). Additional recent balloon measurements corroborate this result (Engel et al., 2017). A recent, improved analysis of the same balloon observations confirms the insignificant (weakly positive) mean age trend above about 24 km, but shows a negative mean age trend below this level (Ray et al., 2014). Indications of a negative mean age trend in the lowermost stratosphere around the year 2000 have also been found from aircraft observations of mean age (Bönisch et al., 2011). Hence, in particular in the stratosphere above about 24 km, climate models results are not consistent with observations regarding trends in the BDC.

The strengthening BDC trend in models is clearly related to greenhouse gas-induced tropospheric warming. This warming strengthens the subtropical jets, thereby shifting the critical levels for wave breaking upwards and equatorwards and intensifying the mechanical forcing of the BDC (e.g. Garcia and Randel, 2008; Shepherd and McLandress, 2011; Garny et al., 2011). A recent study by Oberländer et al. (2016) indicates that the BDC trend is consistent with the expansion of the troposphere and an associated upward shift of the tropopause with climate change.

Concerning BDC changes over decadal periods, horizontal circulation shifts in the latitudinal direction also cause important effects. This has recently been shown for a southward circulation shift from 2002 to 2012 which likely caused a hemispheric dipole mean age change pattern (the age increased in the NH and decreased in the SH) as observed by the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) satellite instrument (Stiller et al., 2017). Ensembles of climate model simulations may include mean age change patterns which are more complex than a global decrease and even have some resemblance to MIPAS observations (Garfinkel et al., 2017). Hence, to simulate observed past BDC changes it is important to correctly represent variability on inter-annual to decadal timescales in the models. For interpreting mean age trends, the consideration of atmospheric mixing processes turns out to be particularly important (e.g. Ray et al., 2010; Ploeger et al., 2015a; Dietmüller et al., 2017).

Meteorological reanalyses combine a global weather forecast model with atmospheric observations through a data assimilation system to provide an optimal estimate of the atmospheric state from the past to present. These reanalyses are based on an unchanged forecast model version and assimilation system to minimise artificial changes in the state variables due to changes in the reanalysis system. However, as the observational data sets included in the assimilation change over time, abrupt state changes may still occur to some degree, rendering reanalysis-based trend studies challenging. Considering more than one reanalysis in such studies increases the reliability of results considerably. The Stratosphere-troposphere Processes And their Role in Climate (SPARC) Reanalysis Intercomparison Project (S-RIP) aims at an inter-comparison of the current generation reanalysis products, as described by Fujiwara et al. (2017). The present paper contributes to this project by comparing the representation of the stratospheric BDC in the three most modern reanalysis products: (i) ERA-Interim from the European Centre for Medium-Range Weather Forecasts (ECMWF), (ii) JRA-55 from the Japanese Meteorological Agency and (iii) MERRA-2 from the National Aeronautics and Space Administration (NASA).

Past analyses have shown that modern reanalyses provide an improved representation of the BDC (e.g. Diallo et al., 2012; Monge-Sanz et al., 2012) compared with older reanalysis products like the ECMWF's ERA-40 reanalysis (e.g. Monge-Sanz et al., 2007). In particular, ERA-Interim has been shown to combine a negative mean age trend throughout most regions of the stratosphere with a weakly positive age trend in the NH above about 24 km, similar to existing balloon observations (Ploeger et al., 2015a). Moreover, the hemispheric dipole age trend pattern as observed by MIPAS from 2002 to 2012 turns out to be reproduced, at least qualitatively, by ERA-Interim-driven simulations (Ploeger et al., 2015b). However, the robustness of these results concerning the representation of the BDC in different reanalyses is an open question.

A very recent study by Chabrillat et al. (2018) compares the BDC in various reanalyses by using a kinematic transport model (Belgian Assimilation System for Chemical ObsErvations, BASCOE) within the scope of the S-RIP project. The results presented here are based on a diabatic transport model (for further details see Sect. 2) and therefore complement the study by Chabrillat et al. (2018) regarding the representation of vertical transport. The main goal of our paper is to assess the robustness of the climatology and seasonality of the BDC as well as its trends in current generation reanalyses, as imprinted on stratospheric age of air. For that reason, we calculate and analyse the mean age of air as well as the full time-dependent stratospheric age spectrum, which has not been used for a model inter-comparison of the BDC to date.

The modelling method (transport model and age of air diagnostics) is described in Sect. 2. In Sect. 3 the climatology and seasonality of age of air from different reanalyses is compared before considering trends in Sect. 4. Section 5 provides a comparison to existing observational mean age estimates, and Sect. 6 discusses the results, in particular compared to the results from the complementary model study of Chabrillat et al. (2018).

2 Data and method

The Chemical Lagrangian Model of the Stratosphere (CLaMS) is a Lagrangian model for calculating transport and chemistry for trace gas species based on the motion of 3-D forward trajectories and an additional parameterised representation of atmospheric small-scale mixing processes (McKenna et al., 2002). Parameterised mixing in the model is driven by deformations in the large-scale flow, such that, in regions of large flow deformations, strong mixing occurs (Konopka et al., 2004). Model transport is calculated in an isentropic vertical coordinate framework with potential temperature θ being the vertical coordinate throughout the stratosphere and upper troposphere, and with the cross-isentropic vertical velocity deduced from the total diabatic heating rate (of the respective reanalysis forecast), including effects of radiative and turbulent heating, as well as latent heat release. Further details about the CLaMS model set-up used in this study can be found in Pommrich et al. (2014).

As described by Ploeger and Birner (2016), and briefly reviewed in the following, we calculate the age of air spectrum for each reanalysis from multiple tracer pulses and the mean age from the spectrum. The age spectrum G is a boundary value Greens function for the continuity equation of a conserved and passive trace gas species (e.g. Hall and Plumb, 1994; Holzer and Hall, 2000) and relates the trace gas mixing ratio χ(r,t) at location r and time t to the mixing ratio χ₀(t) at the boundary surface, where it is assumed to be uniform (e.g. Waugh and Hall, 2002):

$$\begin{array}{}\text{(1)}& \mathit{\chi }(r,t)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\mathit{\tau }{\mathit{\chi }}_{\mathrm{0}}(t-\mathit{\tau })G(r,t,\mathit{\tau }).\end{array}$$

Here, the integration is taken along transit time τ, and the boundary surface is usually taken to be the tropical tropopause or the tropical surface, a particularly common choice in models. The age spectrum calculation applied in this study is set up analogously to the approach described by Ploeger and Birner (2016), and similarly to the calculation in the GEOS climate model (Li et al., 2012a, b). The calculation method is based on N=60 inert pulse tracers, approximating a delta distribution lower boundary condition ${\mathit{\chi }}_{\mathrm{0}}^i\left(t\right)=\mathit{\delta }(t-t_i)$ with $i=\mathrm{1},\mathrm{\dots },\mathrm{60}$, defining tracer pulses at source times t_i. For such a set of pulse tracers, the age spectrum at transit time τ_i is related to the tracer mixing ratio of the ith species via

$$\begin{array}{}\text{(2)}& G(r,t,{\mathit{\tau }}_i)={\mathit{\chi }}_i(r,t);\end{array}$$

hence the age spectrum can be directly calculated from the pulse tracer mixing ratios in the simulation.

To approximate the delta distribution characteristics, pulse tracer mixing ratios are set to one in the lowest (orography following) model layer (approximately the boundary layer) in the tropics between 15^∘ S and 15^∘ N for 30 days. The 60 different species are pulsed every other month (first species in the first month, second species in the third month, etc.), such that after 10 years of simulation all species have been pulsed. Hence, after every 10 years the first species is reset to zero and pulsed again. This boundary impulse (time-)evolving response (BIER) method (Ploeger and Birner, 2016) resolves the age spectrum along 10 years of the transit time axis with a bin size of 2 months.

The mean age of air is then calculated as the first moment of the age spectrum:

$$\begin{array}{}\text{(3)}& \mathrm{\Gamma }(r,t)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\mathit{\tau }\mathit{\tau }G(r,t,\mathit{\tau }).\end{array}$$

As the mean age strongly depends on the age spectrum tail, which generally shows an exponential decay after about 4–5 years (e.g. Li et al., 2012a), the effect of the finite age spectrum tail on the mean age may be corrected by fitting an exponentially decaying function (e.g. Diallo et al., 2012; Ploeger and Birner, 2016). Hence, a corrected age spectrum may be defined by extrapolating the spectrum tail for transit times $\mathit{\tau }>{\mathit{\tau }}^{*}=\mathrm{10}$ using

$$\begin{array}{}\text{(4)}& G_{\mathrm{corr}}(r,t,\mathit{\tau })=G(r,t,{\mathit{\tau }}^{*})e^{-\frac{(\mathit{\tau }-{\mathit{\tau }}^{*})}{\mathit{\xi }(r,t)}},\end{array}$$

with $G(r,t,{\mathit{\tau }}^{*})$ the age spectrum value at 10 years, and the tail decay timescale ξ(r, t) (depending on location and time) being estimated from the exponential fit to the spectrum at transit times between 5 and 10 years. As the full age spectrum is the probability distribution of transit times, it is generally normalised to unity. However, due to the truncation of the simulated spectrum at 10 years, the integration of the spectrum over transit time leads to a norm less than one (see Fig. 6):

$$\begin{array}{}\text{(5)}& N^{\mathrm{10}}=\underset{\mathrm{0}}{\overset{\mathrm{10}\mathrm{year}}{\int }}\mathrm{d}\mathit{\tau }G(r,t,\mathit{\tau })<\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\mathit{\tau }G(r,t,\mathit{\tau })=\mathrm{1}.\end{array}$$

Including the exponential tail correction improves the normalisation, but small differences to unity remain due to the finite resolution along the transit time axis. On the one hand, including the correction generally improves comparisons of mean age to observations. On the other hand, MERRA-2 age spectra show a much more pronounced tail without a clear exponential decay over 10 years in some cases (see Sect. 5), violating the necessary assumption for the finite tail correction. Hence, for most parts of the analysis we simply consider the finite tail age spectra over 10 years without including the tail correction. This simplifies the interpretation of the comparison of different reanalyses by only considering the resolved part of the age spectrum. Effects of the unresolved tail and the finite tail correction are further discussed in Sect. 5.

Two other age spectrum-based transport diagnostics are the modal age, which is the transit time of the maximum age spectrum peak, and the age spectrum width

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }(r,t)=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathrm{d}\mathit{\tau }[\mathit{\tau }-\mathrm{\Gamma }(r,t\left){]}^{\mathrm{2}}G\right(r,t,\mathit{\tau })}.\end{array}$$

The spectrum width is strongly influenced by long transit times in the spectrum tail and therefore is usually considered a measure of the strength of recirculation (e.g. Li et al., 2012a). The modal age, in comparison, can be interpreted as a measure of the residual circulation in the tropics and in the winter stratosphere as it is closely related to the residual circulation transit time in these regions. As additional diagnostics for the interpretation of processes affecting zonal mean mean age we consider residual circulation transit times (RCTTs, Birner and Bönisch, 2011) and the net mixing effect on mean age (“ageing by mixing”, Garny et al., 2014) in Sect. 6. The RCTT is the transit time of a (hypothetical) air parcel if it was transported by the residual circulation alone and, by definition, solely includes effects of the residual circulation (Birner and Bönisch, 2011). For this paper, RCTTs are calculated with the CLaMS trajectory module and using the zonal mean diabatic residual circulation in isentropic coordinates $({\stackrel{\mathrm{‾}}{v}}^{*},{\stackrel{\mathrm{‾}}{Q}}^{*})$ (Ploeger et al., 2015b). Here, Q denotes the diabatic heating rate, v denotes the meridional velocity component, overlined quantities denote mass-weighted zonal averages and primed quantities represent the respective fluctuations therefrom (e.g. Andrews et al., 1987, chap. 9.4); hence, ${\stackrel{\mathrm{‾}}{v}}^{*}=\left(\stackrel{\mathrm{‾}}{\mathit{\sigma }v}\right)/\stackrel{\mathrm{‾}}{\mathit{\sigma }}$, and ${\stackrel{\mathrm{‾}}{Q}}^{*}=\left(\stackrel{\mathrm{‾}}{\mathit{\sigma }Q}\right)/\stackrel{\mathrm{‾}}{\mathit{\sigma }}$.

In addition, the zonal mean mean age $\stackrel{\mathrm{‾}}{\mathrm{\Gamma }}$ is affected by mixing processes in the atmosphere. The local effect of this mixing, the local eddy mixing tendency, is represented in the zonal mean (isentropic) tracer continuity equation for mean age by the divergence ℳ of a 2-D mixing flux vector (e.g. Andrews et al., 1987, Eq. 9.4.21). The net mixing effect on the zonal mean mean age, the ageing by mixing, is then calculated by integrating the local eddy mixing tendency ℳ along residual circulation trajectories x(t) (Ploeger et al., 2015b):

$$\begin{array}{}\text{(7)}& \stackrel{\mathrm{‾}}{\mathrm{\Gamma }}(x,t)={\mathit{\tau }}_{\mathrm{RCTT}}(x,t)+\underset{t_{\mathrm{0}}}{\overset{t}{\int }}\mathcal{M}\left(x\right(t^{\prime }\left)\right)\mathrm{d}{t}^{\prime }.\end{array}$$

Here, τ_RCTT is the RCTT, for transport from the 340 K isentropic surface in the tropics (30^∘ S–30^∘ N), and t₀ is the time when the trajectory intersected the 340 K surface. As proposed by Garny et al. (2014) and recently further evidenced by Dietmüller et al. (2017), the net eddy mixing effect can be well approximated by the difference between the mean age and the RCTT. We apply this approximation in the following without explicitly calculating the mixing effect (see Sect. 6). For further details about the RCTT and the ageing by mixing calculation see Ploeger et al. (2015b).

The CLaMS simulations are driven with horizontal winds and diabatic heating rates from the three most recent reanalysis data sets: ERA-Interim, JRA-55 and MERRA-2. The simulations for ERA-Interim and JRA-55 both start on 1 January 1979, whereas the MERRA-2 simulation starts on 1 January 1980. Due to the specific pulse tracer set-up described above it takes 10 years of simulation until all pulse tracers have been set and the age spectrum can be evaluated. To enable the age of air analysis beginning in 1979 (in 1980 for MERRA-2), for comparison with balloon-borne mean age measurements in Sect. 5, a 10-year-long model spin-up is carried out by repeating conditions of the first simulation year. However, the simulated age of air before 1989 includes the effect of the spin-up and most parts of the analysis presented here are restricted to the time period from 1989 to 2015. For MERRA-2, the year 1989 still includes a very weak spin-up effect, but without any effect on our conclusions, as all results of the paper could be analogously derived for the 1990–2015 period (not shown).

The three different reanalyses used here have been recently described by Fujiwara et al. (2017), with further details given by Dee et al. (2011) for ERA-Interim, by Kobayashi et al. (2015) for JRA-55 and by Gelaro et al. (2017) for MERRA-2. For driving the CLaMS model simulations, reanalysis horizontal winds and diabatic heating rates from the reanalysis forecast are used on native model levels and with a horizontal resolution of $\mathrm{1}{}^{\circ }\times \mathrm{1}$^∘ in latitude and longitude. The age of air results from the different simulations have been interpolated to potential temperature levels (same for all reanalyses) and monthly zonal mean climatologies have been created.

Figure 1 Mean age climatology (1989–2015) for December–February (DJF, a–c) and June–August (JJA, d–f) for ERA-Interim (a, d), JRA-55 (b, e) and MERRA-2 (c, f). Thin solid black lines highlight particular mean age contours, the thin dashed black lines show pressure levels in hPa, and the thick black line is the (lapse rate) tropopause (calculated from each reanalysis following WMO, 1957).

Figure 2 Age spectrum on the 500 K potential temperature surface for December–February in (a) the tropics between 20^∘ S and 20^∘ N, (b) middle latitudes between 30 and 60^∘ N and (c) polar regions between 60 and 90^∘ N. Different colours represent different reanalyses (ERA-Interim is shown using black, JRA-55 is shown using blue and MERRA-2 is shown using red). Vertical dashed lines indicate the respective mean age. (Note the different y axis range for tropical and extratropical spectra.)

3 Seasonal variations in age of air

Climatological mean age data (1989–2015) from the three reanalyses for boreal winter (December–February, DJF) and summer (June–August, JJA) are compared in Fig. 1. Throughout most parts of the stratosphere JRA-55 shows the youngest mean age and MERRA-2 shows the oldest mean age. Values for ERA-Interim lie in between. However, in the tropical lower stratosphere below about 700 K potential temperature (about 20 hPa) ERA-Interim shows the youngest mean age of the three reanalyses, which is consistent with overly strong tropical upwelling as documented in the literature (e.g. Dee et al., 2011).

The MERRA-2 mean age is substantially older than the other reanalyses, with the largest relative differences in the tropical stratosphere (e.g. about 2.5 years for MERRA-2 at about 500 K (about 50 hPa) compared to about 1.5 years for ERA-Interim). Moreover, at high latitudes the MERRA-2 mean age is clearly older, reaching maximum values of more than 5.5 years in the southern polar vortex compared with less than 5 years for ERA-Interim and JRA-55 (note that these values are calculated from the age spectrum over 10 years and that actual mean ages are higher).

Despite these differences in climatological average mean age values, seasonal variations (as estimated from December–February to June–August differences) are very similar in the three reanalyses. The youngest air is found in the tropical stratosphere during boreal winter, when the BDC upwelling maximises. Conversely, the oldest air occurs in the high-latitude stratosphere during wintertime and is related to the strongest downwelling within the deep BDC branch in that region and season. Furthermore, a flushing of the summertime lower stratosphere with young air masses from the tropics (e.g. Hegglin and Shepherd, 2007; Bönisch et al., 2009) is evident in all reanalyses. This summertime flushing shows a robust hemispheric asymmetry, and is stronger in the NH than in the SH (Konopka et al., 2015).

More details about stratospheric transport can be observed in the age spectrum, presented in Fig. 2 for the tropics, middle and high latitudes on the 500 K potential temperature surface. The tropical age spectra show an almost unimodal shape with a clear peak at short transit times, and only very weak additional peaks and a decaying tail at larger transit times. In middle and high latitudes the age spectra show distinct multiple peaks caused mainly by the seasonality of transport into the stratosphere (e.g. Reithmeier et al., 2007; Li et al., 2012a; Diallo et al., 2012; Ploeger and Birner, 2016). These multiple spectrum peaks are robustly found for all reanalyses, again indicating robustness in the representation of seasonal variations in stratospheric transport in modern reanalyses. In general, ERA-Interim and JRA-55 age spectra are very similar for all regions, but MERRA-2 spectra differ more substantially. In particular, for MERRA-2 the secondary peaks at higher ages are delayed by a few months compared with ERA-Interim and JRA-55, again indicating slower transport for MERRA-2. At high latitudes, the modal age (transit time of the maximum spectrum peak) for MERRA-2 may occur with a delay of more than a year (e.g. Fig. 2c). In all regions the modal peak is shifted to larger transit times (higher ages) for MERRA-2. In particular, the spectrum tail is much more pronounced for MERRA-2 (e.g. Fig. 2f), with age spectrum values more than twice as large as ERA-Interim and JRA-55 at transit times larger than about 8 years. Hence, there is a substantially larger fraction of very old air for MERRA-2 than for ERA-Interim and JRA-55.

Figure 3 Age spectrum (as transit time probability density function (PDF) in month⁻¹) at 400 K for December–February (a, b, c) and June–August (d, e, f) for ERA-Interim (a, d), JRA-55 (b, e) and MERRA-2 (c, f). Climatological values for 1989–2015 are shown. The black line shows the mean age, and the white diamonds show the modal age.

Figure 4 Same as Fig. 3, but at the 600 K potential temperature surface.

Figure 5 Cross-isentropic vertical velocity dθ∕dt from the total diabatic heating rate for ERA-Interim (a), JRA-55 (b) and MERRA-2 (c). Climatological annual mean distributions for 1989–2015 are shown. The thin black lines show pressure contours, and the thick black line is the (lapse rate) tropopause.

The global view of the age spectrum in the lowest stratosphere at 400 K is presented in Fig. 3. The 400 K isentrope has been chosen as a representative level for the shallow BDC branch. The robust representation of transport seasonality is again evident from the similar occurrence of the multiple spectrum peaks. In particular, the modal peak (modal age) of the tropical age spectra is shifted to younger ages in boreal winter, which is consistent with faster wintertime BDC upwelling. The strong flushing of the Northern Hemisphere (Bönisch et al., 2009), and to a lesser degree also of the Southern Hemisphere, is clearly visible in the extension of the modal peak (white diamonds) from the tropics into the summer hemisphere in all reanalyses. This extension of the tropical young air signal deep into the summer hemisphere is consistent with the general understanding of a less isolated tropics and stronger isentropic mixing causing horizontal exchange between the tropics and the extratropics during summer.

The subtropics in Fig. 3 are characterised as transition regions between approximately monomodal tropical age spectra and middle- and high-latitude spectra with distinct multiple peaks. The contrast between tropical and extratropical age spectra is strongest in the respective winter hemisphere, which is indicative of a stronger subtropical jet and related transport barrier compared to the summer hemisphere. For MERRA-2, the transition between tropical and extratropical age spectra is more dilute, indicating stronger exchange between tropics and middle latitudes in the lowest stratosphere. This stronger exchange likely results in a stronger recirculation of extratropical aged air masses into the tropics, causing older air throughout the stratosphere.

Very similar results emerge from the comparison of the age spectra at higher levels in the stratosphere (here 600 K in Fig. 4), where the deep BDC branch dominates transport. The seasonality in tropical upwelling (faster in boreal winter) as well as in subtropical transport barriers (stronger in boreal winter) is robustly represented in the different reanalysis data sets. In particular the older modal age for MERRA-2 in the tropics, evident at 400 K (Fig. 3), becomes even more clear at 600 K.

This shift of the tropical modal age to longer transit times is related to a slower tropical upwelling for MERRA-2, which is consistent with the comparison of reanalysis climatological total diabatic heating rates shown in Fig. 5. In particular throughout the tropical tropopause layer (e.g. Fueglistaler et al., 2009a) and up to about 500 K the MERRA-2 heating rates are lower than for the other reanalyses (Fig. 5). Below the tropopause around the level of zero radiative heating (around 360 K) there is even a cooling layer evident in the MERRA-2 heating rates causing a “bottleneck” structure in annual mean upwelling. The differences in heating rates could be related to ozone differences between the reanalyses (e.g. Davis et al., 2017). ERA-Interim uses prescribed climatological ozone fields for the radiative calculations, including a potential high bias in ozone concentrations in the tropical lower stratosphere (e.g. Fueglistaler et al., 2009b). JRA-55, in comparison, uses prescribed time-varying ozone fields from a model simulations (Kobayashi et al., 2015), and MERRA-2 uses interactive ozone for radiation and heating rate calculations (e.g. Wargan et al., 2018).

Regarding the age spectrum, MERRA-2 further shows a weaker gradient between tropical and extratropical age spectra compared with the other reanalyses at 600 K (Fig. 4). This weaker tropical–extratropical contrast indicates stronger exchange between tropics and extratropics and a weaker subtropical transport barrier. Furthermore, the MERRA-2 age spectra at 600 K show a more pronounced spectrum tail, as previously noted. Differences in the age spectrum tail at transit times larger than about 5 years cannot be caused by the differences in the diabatic circulation (in Fig. 5) alone, because related transit times along the residual circulation are generally below about 5 years (e.g. Birner and Bönisch, 2011, Fig. 2). Therefore, the more pronounced spectrum tail for MERRA-2 compared with ERA-Interim and JRA-55 is likely a result of the stronger recirculation of air into the tropics at lower levels. Recirculation into the tropics, particularly at low levels, causes air masses to circulate through the stratosphere within the BDC several times before sinking back into the troposphere and significantly increases the age (e.g. Neu and Plumb, 1999; Garny et al., 2014). Therefore, an increased recirculation enhances the fraction of aged air masses and hence the tail of the age spectrum (e.g. Li et al., 2012a).

Figure 6 further compares the characteristics of the spectrum tail between the three reanalyses by showing the spectrum normalisation (coloured shading) and the tail decrease timescale (blue contours). The spectrum norm was calculated as the integral of the age spectrum over 10 years (Eq. 5), and the tail decrease timescale was calculated as the exponential decrease rate of the spectrum tail at transit times larger than 5 years (Eq. 4). As the model age spectra are truncated after 10 years (see Sect. 2) a stronger spectrum tail will result in a less stringent normalisation condition and in a larger difference of the age spectrum norm from unity. Indeed, Fig. 6 shows that the MERRA-2 age spectra are less well normalised, with a norm below 0.8 throughout large regions of the stratosphere, whereas for ERA-Interim and JRA-55 the norm is always above 0.9 (in the lower part of the stratosphere even above about 0.95). This lack of normalisation of the MERRA-2 age spectra is clearly related to a much larger tail decrease timescale (global average 5.13 years) compared with ERA-Interim (2.84 years) and JRA-55 (2.91 years). Compared to the climate model simulated age spectra of Li et al. (2012a), with a global average tail decay timescale of 2.77 years, ERA-Interim and JRA-55 show comparable values, whereas MERRA-2 differs substantially.

Figure 6 Norm of the climatological (1989–2015) age spectra for ERA-Interim (a), JRA-55 (b) and MERRA-2 (c). Blue contours show the climatological tail decay timescale (in years), the thin black dashed lines show the pressure levels in hPa, and the thick black line is the (lapse rate) tropopause.

Figure 7 Long-term mean age trends for 1989–2015 (a, b, c), decadal trends for 1989–2001 (d, e, f) and decadal trends for 2002–2015 (g, h, i), for ERA-Interim (a, d, g), JRA-55 (b, e, h) and MERRA-2 (c, f, i). The trend significance, estimated in multiples of the standard deviation σ, is shown as grey contours (2σ contour thick grey, and then decreasing in 0.2 steps as thin lines). The thin black solid contours show the climatological mean age distribution, the thin black dashed lines show pressure levels in hPa, and the thick black line is the (lapse rate) tropopause.

4 Age of air trends

In the following, we analyse trends in age of air for the same periods as considered by Chabrillat et al. (2018) to simplify the comparison to their results in Sect. 6. We refer to 1989–2015 as “long-term” trends and to 1989–2001 (“pre-2000” in the following) and 2002–2015 (“post-2000”) as “decadal” trends. It should be noted that even the long-term period spans only 27 years, which is relatively short compared with climate model simulation periods but is the longest period we can obtain from the reanalyses without including spin-up effects in the results (see Sect. 2). Throughout the paper, trends are calculated from the linear regression of monthly mean time series after deseasonalising (by subtracting the mean annual cycle), and the significance is measured in multiples of the standard deviation of the linear trend.

4.1 Long-term trends from 1989 to 2015

Long-term mean age trends for 1989–2015 are presented in Fig. 7a–c. Overall, mean age trends are negative throughout most regions of the lower stratosphere below about 600 K (30 hPa) and for all three reanalyses. Hence, the representation of long-term changes in mean age in the lower stratosphere in current reanalyses is largely robust and indicates an accelerating shallow BDC branch. Mean age trends above about 600 K are not robust, with ERA-Interim showing positive trends in the NH whereas JRA-55 and MERRA-2 show negative trends. Hence, trends in the deep BDC branch appear not to be robustly represented in current reanalysis data sets. The mean age decrease in the shallow BDC branch is in agreement with an acceleration of the residual circulation in reanalyses (Abalos et al., 2015; Miyazaki et al., 2016).

Closer examination of the 1989–2015 trends in Fig. 7a–c reveals detailed differences between the reanalyses. First, ERA-Interim shows an inhomogeneous pattern with the strongest negative trends in the southern subtropics, and even weakly positive trends in the NH above about 30 hPa. Second, JRA-55 provides a very homogeneous picture, showing negative mean age trends everywhere. Third, MERRA-2 again shows largely negative trends, maximising in the southern subtropics similarly to ERA-Interim, and positive trends in the NH lowermost stratosphere below about 400 K.

Figure 8 Age spectrum trend (annual mean) at 400 K for ERA-Interim (a, d, g), JRA-55 (b, e, h) and MERRA-2 (c, f, i), and for the periods 1989–2015 (a, b, c), 1989–2001 (d, e, f), and 2002–2015 (g, h, i). The thin black contours show the climatological age spectrum, the thick black line shows the climatological mean age, and the white diamonds show the climatological modal age.

Changes in the full age spectrum can shed more light on the processes involved during each period. Figure 8a–c shows age spectrum changes (coloured shading) during the 1989–2015 period in the lower stratosphere at 400 K and at all latitudes for the different reanalyses, overlaid with contours of climatological mean age (thick black lines). Remarkably, trends in the age spectra at young transit times are very consistent, particularly in the tropics, with all reanalyses showing an increase in the fraction of air masses with transit times younger than the modal age, and a simultaneous decrease in the air mass fraction with transit times just older than modal age. These changes indicate a shift of the spectrum peak (modal age) to younger transit times over time. Modal age can be related to the residual circulation transit time (e.g. Li et al., 2012a), particularly in the tropics and the winter hemisphere stratosphere (Ploeger and Birner, 2016). The shift of the modal age to younger ages in the tropics indicates an acceleration of the residual mean mass circulation.

This decrease in the modal age emerges even more clearly at higher levels (e.g. 600 K in Fig. 9a–c). Independent of reanalysis and latitude, the spectrum peak shifts to younger transit times over the 1989–2015 period. Differences in mean age trends between the different reanalyses appear to be related to differences in the age spectrum tail. Clearly, the weakly increasing mean age in the NH above about 30 hPa (Fig. 7) in ERA-Interim is related to an increasing age spectrum tail at transit times older than about 4 years (Fig. 9a), which is absent in the other two reanalyses. Note that the 600 K potential temperature level is at the lower boundary of the region of increasing mean age in ERA-Interim, and that the described increase of the age spectrum tail becomes even clearer at levels above. However, 600 K is chosen as it is consistent with the later 2002–2015 trends, which show the mean age dipole change pattern only below about 20 hPa.

In agreement with the general shift of the age spectrum peak towards younger transit times from 1989 to 2015 the fraction of young air masses, younger than 6 months, robustly increases in the lower stratosphere in all three reanalyses (Fig. 10). Differences occur in the NH lowermost stratosphere, where the young air mass increase is strongest for JRA-55 and is absent for MERRA-2, which is consistent with the respective changes in mean age (Fig. 7). Note that the young air mass fraction is a particularly robust diagnostic as it is independent of the truncation of the age spectrum to 10 years. Observable changes in the young air mass fraction are confined below about the 500 K potential temperature level. Above this level, the air is generally older than 6 months, and the young air mass fraction with transit times shorter than 6 months almost vanishes. Changes in the old air mass fraction, older than 2 years, show negative changes in the tropical stratosphere above about 450 K (not shown), which is consistent with the shift of the age spectrum towards shorter transit times (also at higher levels).

Figure 9 Same as Fig. 8, but at the 600 K potential temperature surface.

Figure 10 Trend in the young air mass fraction (transit time less than 6 months) from 1989 to 2015 for (a) ERA-Interim, (b) JRA-55 and (c) MERRA-2, in percent per decade. The white contours show the young air mass fraction climatology in percent, and the thick white line is the tropopause.

4.2 Decadal changes during 1989–2001 and 2002–2015

The decadal mean age changes over shorter periods are much more diverse than the long-term trends. Furthermore, certain characteristics in these decadal changes strongly depend on the start and end points of the period considered and should not be taken as representative of long-term trends. Nevertheless, as observational data sets only exist for restricted periods and several past studies have focused on circulation changes before and after the year 2000, we include a discussion of such decadal changes in the following. The decadal periods considered (1989–2001 and 2002–2015) were chosen for better comparability with Chabrillat et al. (2018).

For 1989–2001, ERA-Interim and MERRA-2 show mean age trend patterns similar to the longer 1989–2015 period (Fig. 7d and f). JRA-55, in contrast, shows much stronger negative trends, particularly in the tropics (Fig. 7e). These strong negative trends in JRA-55 are related to a stronger effect of the volcanic eruption of Mt. Pinatubo in June 1991, which increased mean age in the reanalysis in the lower stratosphere during the 2–3 years following the eruption (Diallo et al., 2017). As the Pinatubo eruption and related ageing of air occurs close to the beginning of the trend period, 1989–2001 negative decadal mean age trends appear enhanced for JRA-55 (see Fig. 11b). More generally, the strong negative mean age trends in the Southern Hemisphere (SH) during the 1989–2001 period are likely related to the effects of ozone depleting substances (Polvani et al., 2018), with the dominant effect caused by their chemical impact with respect to depleting ozone (Abalos et al., 2019). Reanalyses may include ozone depletion effects by assimilating observed temperatures even without having realistic ozone fields.

Mean age trends for the later period 2002–2015 are least consistent among the three reanalyses (Fig. 7g–i). MERRA-2 shows negative trends similar to the 1989–2015 and 1989–2001 periods. Conversely, JRA-55 shows positive trends throughout almost the entire stratosphere, with the exception of a small region in the NH lower stratosphere of insignificant trends. In contrast, ERA-Interim shows decreasing age in the lowest stratosphere, and a clear dipole pattern above with increasing age in the NH lower stratosphere, and decreasing age in the SH. This dipole pattern has been shown to be consistent with mean age trends based on satellite observations (e.g. Haenel et al., 2015; Ploeger et al., 2015b) and observed HCl trends (Mahieu et al., 2014), although for a slightly different period (2002–2012). Recently, Stiller et al. (2017) related the dipole trend pattern to a southward shift of the global stratospheric circulation. Obviously, only ERA-Interim shows this dipole pattern in mean age trends for 2002–2015, whereas the other two reanalyses provide only very weak indications for slightly more ageing in the NH (JRA-55) or a weaker age decrease in that region (MERRA-2).

The differences between the reanalyses are larger for the age spectrum changes over the shorter periods, 1989–2001 and 2002–2015, than for the long-term, 1989–2015, trend. At the lower potential temperature level of 400 K the spectrum peak may shift towards younger transit times (e.g. JRA-55 for 1989–2001) or towards older transit times (e.g. MERRA-2 for 2002–2015), and the spectrum tail may decrease or increase depending on the reanalysis (Fig. 8d–i). At the upper level of 600 K, more consistent spectrum changes emerge. During the early period 1989–2001 all reanalyses show a shift of the spectrum peak towards younger transit times (Fig. 9d–f), which is consistent with the decreasing mean age over this period. During the later 2002–2015 period MERRA-2 mainly shows a shift of the spectrum peak to shorter transit times (Fig. 9i), which is consistent with decreasing mean age (Fig. 7i). JRA-55, in contrast, shows a shift of the spectrum peak towards longer transit times and an increase in the spectrum tail (Fig. 9h), resulting in positive mean age trends (Fig. 7h). ERA-Interim shows a hemispheric asymmetric pattern (Fig. 9g). In the NH, the spectrum peak shifts towards longer transit times and the spectrum tail increases, resulting in an increasing mean age during this period (Fig. 7g). In the SH, opposite changes occur, with a spectrum peak shift towards shorter transit times and a decrease of the spectrum tail. Hence, the hemispheric asymmetric age spectrum changes in ERA-Interim for transit times less than 10 years are consistent with the hemispheric dipole trend in mean age. Clearly, spectrum changes in the tail at transit times of around 5 years appear critical for the positive mean age trend in the NH. As 5 years is beyond the timescale of the residual circulation, mixing effects are likely involved in causing the hemispheric dipole pattern of mean age changes, as recently concluded from analysis of different diagnostics (Ploeger et al., 2015b; Stiller et al., 2017). Interestingly, age spectrum changes in the NH from 2002 to 2015 appear qualitatively consistent between ERA-Interim and JRA-55, showing decreasing values at transit times less than about 2 years and increasing values at transit times around 5 years. Hence, the more pronounced NH ageing in ERA-Interim compared with JRA-55 is the result of a very subtle balance between changes in the age spectrum at short transit times (less than about 2 years) and in the tail (around 5 years).

5 Comparison to observations

To assess the reliability of the representation of the stratospheric BDC in different reanalyses, the mean age of air is compared to the mean age estimated from observations of the long-lived trace gas species SF₆ and CO₂. Figure 11a compares latitude sections of mean age from reanalysis, airborne in situ observations (same data as shown in Waugh and Hall, 2002, based on various measurements; Boering et al., 1996; Andrews et al., 2001; Elkins et al., 1996; Ray et al., 1999; Harnisch et al., 1996) and MIPAS satellite observations of SF₆ (Stiller et al., 2012; Haenel et al., 2015) at 500 K potential temperature (20 km for the observations).

Overall, the reanalysis mean age lies within the uncertainty range of the observations. In the tropics, the in situ observed mean age is significantly lower than the reanalysis age, with only ERA-Interim reaching similar values. However, MIPAS shows a much older mean age (above 2 years), which is similar to MERRA-2 values. Compared to the in situ observations, age gradients in the subtropics are too weak for all reanalyses due to the higher tropical age values. In the middle and high latitudes, ERA-Interim and JRA-55 agree closely with CO₂-based mean age observations, whereas MERRA-2 agrees better with the SF₆-based mean age.

Figure 11 (a) Latitude section of mean age at 20 km from in situ observations (e.g. Waugh, 2009) shown using black symbols, from different reanalysis data sets at 500 K potential temperature (red represents ERA-Interim, blue represents JRA-55 and green represents MERRA-2) and from MIPAS satellite observations (grey). Dashed lines show the reanalysis mean age including the correction for the finite age spectrum tail (see text). Grey shading shows the range between maximum and minimum MIPAS observations at each latitude. (b) Mean age time series in the northern middle latitudes (40–50^∘ N and 30–5 hPa, 600–1200 K for reanalyses). Coloured lines show the mean age from the different reanalyses (red dashed line denoted the corrected mean age for ERA-Interim). Black symbols show the mean age from the balloon observations of Engel et al. (2017), with error bars representing the uncertainty of the observations.

As the reanalysis age spectrum is truncated at 10 years, the respective mean age is biased low by definition. This low bias can be corrected by applying a correction for the finite age spectrum tail (see Sect. 2), and the corrected mean ages are also compared to the observations in Fig. 11a (dashed lines). For all reanalyses the tail correction increases the mean age. This effect is moderate for ERA-Interim and JRA-55, such that the corrected reanalysis mean age remains within the observational uncertainty range. For MERRA-2, in contrast, the effect of the tail correction is large, increasing the mean age by more than 2 years at high latitudes. As a consequence, the tail-corrected MERRA-2 mean age is clearly out of the observational uncertainty range (the tail correction effects will be discussed further below).

Figure 12 (a) Age spectrum at 46^∘ N and 700 K potential temperature for June 1992 from ERA-Interim (a) and MERRA-2 (b). Red lines illustrate the correction method for the finite age spectrum tail by fitting an exponential function to the tail at transit times older than 5 years (see Sect. 2). The dashed red vertical line shows the mean age for the tail-corrected spectrum, the dashed black line shows the mean age for the uncorrected case. (Note the logarithmic y axis, varying over 3 orders of magnitude).

A long-term observational mean age time series only exists for NH middle latitudes from balloon-borne measurements of SF₆ and CO₂ (Engel et al., 2009, 2017). Long-term time series of reanalysis mean age in the NH middle latitudes (averaged between 40 and 50^∘ N and 600 and 1200 K, which is approximately equivalent to the vertical layer for the observations of 30–5 hPa) are compared against this balloon observation time series in Fig. 11b (note, that data points before 1989 include spin-up effects, see Sect. 2). Clearly, the uncertainty arising from using different reanalyses is of a similar magnitude as the uncertainty in the observations, such that no conclusion is possible regarding which reanalysis scenario is most realistic. However, The mean age trends substantially differ between the reanalyses with negative trends for JRA-55 ($-\mathrm{0.17}\pm \mathrm{0.01}$ yr decade⁻¹) and MERRA-2 ($-\mathrm{0.18}\pm \mathrm{0.01}$ yr decade⁻¹) and a positive trend for ERA-Interim (0.07±0.01 yr decade⁻¹), Hence, ERA-Interim appears to agree best with the observed (non-significant) trend of 0.15±0.18 yr decade⁻¹ by Engel et al. (2017). A “best fit” to the observational time series results from applying the finite tail correction to the ERA-Interim mean age, which is in remarkably good agreement with the observations, regarding absolute values, variability and trend (red dashed line).

Figure 11a shows a much stronger sensitivity to the tail correction for the MERRA-2 mean age compared with the other reanalyses. This stronger sensitivity for MERRA-2 results from differences in the spectrum tail compared with the other two reanalyses (see Fig. 2), with a substantially slower tail decrease for MERRA-2. Figure 12 illustrates a case of extreme differences between the reanalyses (June 1992 at 700 K potential temperature and 46^∘ N). For ERA-Interim, the age spectrum values decrease substantially over the 10 years of transit time (more than an order of magnitude), indicating a strongly decaying spectrum tail. For MERRA-2, in contrast, the tail decrease is much slower, leading to a substantially larger tail decay timescale and to a less strict normalisation, as previously discussed above (Fig. 6). This longer tail decay timescale causes the finite tail correction applied to the age spectra resolved over 10 years to have a much larger effect for MERRA-2 than for ERA-Interim and JRA-55. Figure 12 further shows that for MERRA-2 even the assumption of an exponentially decaying age spectrum tail after about 5 years, which has generally been found to be valid for observations (Ehhalt et al., 2004) and models (Li et al., 2012a), is violated in some cases. Hence, the exponential finite tail correction is not applicable for the MERRA-2 age spectra over 10 years of transit time. For that reason, the inter-comparison presented in this paper is based on the uncorrected age spectrum.

6 Discussion

Potential long-term changes in the stratospheric Brewer–Dobson circulation and their relation to anthropogenically forced climate change have been the subject of intense debate over recent years. Climate models show a significant strengthening and acceleration of the BDC (Butchart, 2014), resulting in a global negative mean age trend (e.g. Butchart et al., 2010). Trace gas observations, in comparison, show no indication of a significant long-term mean age trend (Engel et al., 2009, 2017). Studies based on reanalysis meteorology have also provided a diverse picture of mean age trends (e.g. Diallo et al., 2012; Monge-Sanz et al., 2012; Ploeger et al., 2015a) which is partly related to the different trend periods considered.

In view of such diversity in past published mean age trends, the qualitative agreement between the “long-term” trends for 1989–2015 in the lower stratosphere (below about 30 hPa) presented here for the newest three reanalyses is remarkable (Fig. 7). In general, the mean age trend is negative throughout the lower stratosphere with the largest values in the tropical and SH lower stratosphere, indicating an accelerating shallow BDC branch. These negative long-term mean age trends in current reanalyses appear to be largely consistent with mean age trends from climate models over even longer periods (century), representing the forced response of the BDC to increasing greenhouse gas levels. Consistent with these negative mean age trends and the related shift of the age spectrum peak (modal age) towards shorter transit times (e.g. Fig. 9), most current reanalyses show an increasing tropical mass flux (Abalos et al., 2015) and decreasing residual circulation transit times (Thomas Birner, personal communication, 2018). However, Abalos et al. (2015) further pointed out that residual circulation trends may strongly depend on the method used for calculating upwelling velocities.

Above about 30 hPa mean age trends and related circulation changes depend more strongly on the respective reanalysis considered. At these upper levels ERA-Interim shows ageing in the NH, which is not included in the other two reanalyses. Hence, long-term changes in the deep BDC branch appear to be less robustly represented in the reanalyses. The NH ageing in ERA-Interim turns out to be related to a strengthening age spectrum tail indicating changes in recirculation. The evidence for recirculation changes agrees with the results of Ploeger et al. (2015a), which showed that atmospheric mixing processes substantially affect mean age trend patterns, with the strongest effects on shorter timescales such as decades.

Over shorter periods of around a decade the reanalysis age of air trends show substantial differences. For the pre-2000 period, trends in the NH above about 600 K (approximately 24 km) are not consistent (Fig. 7) due to insignificant ERA-Interim trends (Fig. 7d). The strong negative mean age trends in the SH and tropics, in contrast, consistently emerge from all reanalyses and likely represent the combined effects of ozone depleting substances to decrease mean age (Polvani et al., 2018) and of volcanic aerosol from the Pinatubo eruption to increase mean age at the beginning of the 1989–2001 period (Diallo et al., 2017). In particular, regarding the effect of Pinatubo aerosol, differences in the response of the BDC and related differences in age of air between the reanalyses are evident from the middle latitude mean age time series in Fig. 11b. The increase in age after the Pinatubo eruption is very strong in JRA-55, is weaker in ERA-Interim and is absent in MERRA-2. These differences point to different representations of stratospheric volcanic aerosol in the reanalyses causing the differences in decadal BDC trends. In ERA-Interim and JRA-55 the effects of stratospheric volcanic aerosol are only included by the assimilation of observed temperature and wind data, as discussed in more detail by Diallo et al. (2017), whereas MERRA-2 additionally assimilates aerosol optical depth (Fujiwara et al., 2017). Furthermore, the effect of ozone depleting substances (ODSs) is very likely differently represented in the reanalyses due to the usage of different ozone products. ERA-Interim uses an ozone climatology, JRA-55 uses time-varying ozone fields from another model simulation and MERRA-2 uses interactive ozone (see Sect. 3). Hence, differences in the related effects on the BDC are to be expected.

Figure 13 Mean age from ERA-Interim (a), and differences to JRA-55 (b) and MERRA-2 (c), for the 1989–2015 climatological mean. The same for residual circulation transit time RCTT is shown in (d)–(f), and for ageing by mixing in (g–i). The thin black solid contours show the climatology values for the respective reanalysis (e.g. ERA-Interim mean age in a, MERRA-2 mean age in c), the thin black dashed lines show pressure levels, the thick black solid line show the (WMO) tropopause and the blue contours in (e)–(f) show relative differences for RCTTs in percent (solid for positive and dashed for negative differences).

For the post-2000 period, the age of air trends are even less consistent and no clear common circulation change emerges from the comparison of the different reanalyses (Fig. 7). Remarkably, the largest differences between the reanalysis trends occur for the most recent period. This is against the general expectation that a larger consistency between reanalyses emerges for more recent periods, when observational data and the assimilation procedures are more consistent. The inconsistency in the post-2000 trends is mainly related to the fact that JRA-55 and MERRA-2 show opposite mean age trends during the 2002–2015 period throughout the stratosphere (Fig. 7h, i), and that only ERA-Interim shows a hemispheric dipole pattern (Fig. 7g). The hemispheric dipole pattern in the mean age changes in ERA-Interim is related to changes in the age spectrum tail (increasing tail in the NH and decreasing tail in the SH), pointing to changes in recirculation that are in agreement with Ploeger et al. (2015a). As the hemispheric dipole age trend pattern has been recently related to a southward circulation shift by about 5^∘ (Stiller et al., 2017), post-2000 age trend differences between the reanalyses could be caused by differences in the representation of this circulation shift.

Comparison to the recently published reanalysis mean age trends by Chabrillat et al. (2018) may reveal interesting differences regarding the representation of the BDC between kinematic transport models as used in their study (BASCOE model) and diabatic transport models as used here (CLaMS). Regarding the climatology, the mean age from MERRA-2 shows substantially higher values than the other reanalyses for both diabatic and kinematic transport (see Fig. 1, and Chabrillat et al., 2018, Fig. 3). The differences in diabatic heating rates, with weaker tropical upwelling for MERRA-2 (Fig. 5), are qualitatively consistent with the differences in mean age in the diabatic model calculation. The fact that similar mean age differences were found by Chabrillat et al. (2018) suggests that differences in the reanalysis vertical winds, which drive their kinematic model calculations, could be similar to the differences in heating rates.

However, as previously earlier, the differences in heating rates alone cannot fully explain the simulated mean age differences as shown in the following. Figure 13 shows the ERA-Interim climatological mean age and residual circulation transit times (RCTTs), and differences to the other reanalyses. Here, RCTTs are calculated from the diabatic residual circulation in isentropic coordinates, as described in Sect. 2. The mean age from MERRA-2 is older than ERA-Interim by more than 2 years throughout most parts of the stratosphere (Fig. 13c), whereas differences in RCTTs are below 0.5 years (Fig. 13f). Hence, differences in the representation of mixing processes must play a role.

The net mixing effect on mean age (ageing by mixing), calculated as the difference between the mean age and the RCTT following Garny et al. (2014) and further described in Sect. 2, is shown in Fig. 13g–i. This ageing by mixing is about 2 years larger for MERRA-2 than for ERA-Interim in most parts of the lower stratosphere, and clearly explains the large difference in the mean age between the reanalyses. However, ageing by mixing depends on both local mixing characteristics in the lower stratosphere and the transit time of tropical upwelling, which controls the timescale for mixing to affect the ascending air (see Garny et al., 2014; Ploeger et al., 2015a). Hence, the differences in the ageing by mixing between ERA-Interim and MERRA-2 may be caused by either differences in local mixing (eddy diffusivity) or by differences in the RCTTs. From a tropical leaky pipe model perspective ageing by mixing is even linearly related to RCTT, with a longer RCTT causing a larger ageing by mixing (Neu and Plumb, 1999; Dietmüller et al., 2018). In percent, the differences in RCTTs between MERRA-2 and ERA-Interim in the tropical lower stratosphere are large (about 50 %, see blue contours in Fig. 13f), and likely cause a related difference in ageing by mixing. Overall, the larger ageing by mixing for MERRA-2 compared with ERA-Interim clearly causes the differences in mean age. Deeper insight into potential differences in local mixing characteristics could be gained by a reanalysis inter-comparison of effective diffusivity (e.g. Haynes and Shuckburgh, 2000; Abalos et al., 2016).

For the long-term mean age trends (1989–2015) both diabatic and kinematic transport representations result in very similar results for all three reanalyses, as evident from the comparison of Fig. 7 with Fig. 12 of Chabrillat et al. (2018). This good agreement regarding long-term trends again points to a robust representation of the effect of greenhouse gas-induced warming in forcing trends in the BDC. However, regarding decadal changes over shorter periods, the results from diabatic and kinematic transport models show more substantial differences. The strongest differences emerge for the pre-2000 period for ERA-Interim and for the post-2000 period for MERRA-2. For the former case, diabatic transport in CLaMS results in negative mean age trends throughout the SH and tropical stratosphere (Fig. 7d), whereas kinematic transport in BASCOE results in positive trends in these regions (Chabrillat et al., 2018, Fig. 12). As pre-2000 decadal trends are likely controlled by the effects of ozone depleting substances and Pinatubo volcanic aerosol on the BDC, the differences in age trends for this period point to differences in the representation of these processes between diabatic and kinematic transport. In particular differences in the effect of volcanic aerosol (e.g. due to the Pinatubo eruption in June 1991), with increasing mean age for diabatic transport and no significant age changes for kinematic transport seem to be critical, and have also been noted by Chabrillat et al. (2018). For the latter case (post-2000 period for MERRA-2), diabatic transport in CLaMS results in negative trends throughout the stratosphere (Fig. 7i), whereas kinematic transport in BASCOE results in positive age trends (Chabrillat et al., 2018, Fig. 12). Furthermore, less striking differences between diabatic and kinematic mean age trends emerge for MERRA-2 during the pre-2000 period, where kinematic age trends are more negative compared with diabatic trends (Fig. 7f and Chabrillat et al., 2018, Fig. 12). Weak differences also occur for ERA-Interim from 2002 to 2015, where the NH ageing appears stronger in the diabatic age trends (Fig. 7g and Chabrillat et al., 2018, Fig. 12). The causes for the differences regarding decadal trends in the BDC on shorter periods between diabatic and kinematic transport representations, as well as between the different reanalyses, are unclear to date, but appear to be a promising subject for future analysis aiming to improve consistency between the reanalyses.

7 Conclusions

We compared stratospheric mean age and the full age spectrum from simulations with the diabatic CLaMS model driven by different reanalyses (ERA-Interim, JRA-55 and MERRA-2) to investigate the robustness in the representation of the climatology, seasonality and trends of the Brewer–Dobson circulation in current generation reanalysis data sets. Considering the full (time-dependent) age spectrum in a data inter-comparison is novel and provides clearer insights into circulation differences than considering mean age alone. Climatological mean age differs significantly between the different reanalyses, with JRA-55 showing the youngest age and MERRA-2 showing the oldest age throughout almost the entire stratosphere. The substantially older mean age for MERRA-2 appears to be related to a more pronounced age spectrum tail, indicating a stronger effect of the recirculation of air into the tropics. A comparison of residual circulation and mixing effects on mean age confirms that the net mixing effect (including recirculation) is necessary to explain the mean age differences between MERRA-2 and ERA-Interim. The seasonality in the BDC is robustly represented, with very similar mean age seasonality and similar seasonal age spectrum peaks emerging for all reanalyses. Comparison to balloon-borne mean age observations reveals a similarly large spread in simulated and observed mean ages and allows no clear conclusions to be drawn regarding the reliability of the different reanalyses.

In particular, long-term trends in the lower stratosphere (below about 30 hPa) during the 1989–2015 period are robustly represented in the reanalyses mainly showing decreasing mean age, which is strongest in the SH and tropical lower stratosphere. Related to this mean age decrease is a robust shift of the age spectrum towards shorter transit times and an increase in the fraction of young air masses. These long-term age of air changes from reanalyses resemble results from climate model simulations over even longer periods, which simulate an acceleration of the shallow branch of the stratospheric BDC as the forced response to increasing greenhouse gas levels. At upper levels (above 30 hPa), mean age changes in the reanalyses appear less robust, pointing to a less robust representation of changes in the deep BDC branch.

For shorter periods of about a decade (here 1989–2001 and 2002–2015), the age of air changes are more diverse and depend on the specific reanalysis considered. These decadal age changes may even disagree in sign globally for certain periods. Moreover, the hemispheric asymmetric dipole in mean age trends for 2002–2015, as viewed by satellite observations, only emerges for ERA-Interim. Hence, decadal variability in the Brewer–Dobson circulation and the various factors involved (e.g. QBO, ENSO, ODS and volcanic aerosol) turn out not to be robustly represented in current generation reanalyses.