Age of stratospheric air in the ERA-Interim

. The Brewer-Dobson mean circulation and its variability are investigated in the ERA-Interim over the period 1989-2010 by using an off-line Lagrangian transport model driven by analysed winds and heating rates. At low and mid-latitudes, the mean age of air in the lower stratosphere is in good agreement with ages derived from aircraft, high altitude balloon and satellite observations of long-lived tracers. At high latitude and in the upper stratosphere, we ﬁnd, however that the ERA-Interim ages exhibit an old bias, typically of one to two years. The age spectrum exhibits a long tail except in the low tropical stratosphere which is modulated by the annual cycle of the tropical upwelling. The distribution of ages and its variability is consistent with the existence of two separate branches, shallow and deep, of the Brewer-Dobson circulation. Both branches are modulated by the tropical upwelling and the shallow branch is also modulated by the subtropical barrier.Thevariability of


Introduction
Over the last twenty years, the Brewer-Dobson circulation has been recognized as a major component of the climate system (Andrews et al., 1987;Holton et al., 1995;Callaghan, 2005, 2006) which affects radiative budget and atmospheric circulation.
Reanalysed winds from operational weather centres are used to drive Chemistry Transport Models (CTM). Therefore, they are required to properly represent the Brewer-Dobson circulation in order to account for the dependence of the distribution of chemical species in the stratosphere onto the transport properties. It is also important per se to assess the ability of the combined system of a numerical weather forecast model and the associated assimilation system to reproduce the observed behaviour of the stratospheric circulation.
A commonly used metric of the Brewer-Dobson circulation is the age of air, defined as the time spent by a particle in the stratosphere since its entry across the tropopause (Li and Waugh, 1999;Waugh and Hall, 2002). As each air parcel is a mixture of particles with different histories and ages, the age of the parcel is an average over these particles (Kida, 1983;Hall and Plumb, 1994). The age can be further

Published by Copernicus Publications on behalf of the European Geosciences Union.
averaged over time and space to define a mean age over this ensemble or can be described as a distribution denoted as the age spectrum (Waugh and Hall, 2002). A main advantage of the age of air is that it can be estimated from observations of long-lived species (Andrews et al., 2001b;Waugh and Hall, 2002;Stiller et al., 2008;Garcia et al., 2011). The age of air is also used as a mean to compare models (Eyring et al., 2006). The distribution of mean age in latitude and altitude is convenient to visualize the Brewer-Dobson circulation, its strength and its variability (Li and Waugh, 1999;Waugh and Hall, 2002;Austin and Li, 2006). The age spectrum further characterizes the variability and the distribution of transport paths and mixing in the stratosphere (Andrews et al., 1999;Schoeberl et al., 2003Schoeberl et al., , 2005Reithmeier et al., 2008;Li et al., 2012).
Another metric of the Brewer-Dobson circulation is based on the calculation of the residual vertical and meridional velocities (Andrews et al., 1987) which are a representation of the mean zonally averaged mass transport in the stratosphere. This residual circulation is used to calculate transit times from the tropopause crossing (see, e.g. Birner and Bonisch, 2011). Transit times, however, are generally not identical to the age of air as this latter is also influenced by fast stirring and mixing induced by horizontal quasi-isentropic motion in the stratosphere (Waugh and Hall, 2002;Birner and Bonisch, 2011).
The Brewer-Dobson circulation undergoes an annual cycle and changes from year to year. A major mode of variability is the quasi-biennial oscillation (QBO) (Baldwin et al., 2001) which triggers a modulation of vertical transport in the stratosphere by affecting temperature and thus heating rates (Niwano et al., 2003;Punge et al., 2009). An other important factor are volcanic eruptions: in 1991, the Pinatubo has injected massive amount of dust in the stratosphere which have affected its circulation for several years (Thompson and Solomon, 2009). ENSO (Shu et al., 2011) and solar variations are two other sources of Brewer-Dobson variability.
A major source of concern is the existence of a trend in the Brewer-Dobson circulation. Changes in wave propagation and dissipation, as well as possible increases in tropospheric wave activity, are thought to be the primary driver of a strengthened Brewer-Dobson circulation obtained in many models (Butchart and Scaife, 2001;Sigmond et al., 2004;Butchart et al., 2006;Li et al., 2008;Garcia and Randel, 2008). Thompson and Solomon (2005) have observed a cooling of the tropical stratosphere in radiosonde records over the last decades, which is consistent with increased upwelling in the tropical stratosphere.
The analysis of tracer data, however, does not provide evidence for such trend. Engel et al. (2009) and Stiller et al. (2012) even suggest than the age of air might be increasing in some parts of the stratosphere. A possible reason for this discrepancy is that short term trends over one decade are not representative of the trend over one century (Waugh, 2009). Another possibility, however, is that the models or the diag-nostics do not fully account for the stratospheric processes. Ray et al. (2010) analysed the observed trends in mean age and ozone, assuming a simple tropical pipe model, and concluded that "the best quantitative agreement with the observed mean age and ozone trends over the past three decades is found assuming a small strengthening of the mean circulation in the lower stratosphere, a moderate weakening of the mean circulation in the middle and upper stratosphere, and a moderate increase in the horizontal mixing into the tropics". Similarly, Bonisch et al. (2011) found an increase of the Brewer-Dobson circulation in the lower stratosphere but no change at upper levels. In another recent study, Monge-Sanz et al. (2012) found a small old trend of the ages in the same ERA-Interim data used in the present study.
In this study, we present the age of stratospheric air over the period 1989-2010 from Langrangian transport calculations based on most recent reanalysed winds and heating rates from the ERA-Interim reanalysis of the European Centre for Medium Range Weather Forecast (ECMWF). The age of stratospheric air is investigated using backward deterministic trajectories which are integrated over 10 yr in time to evaluate the residence time in the stratosphere. We describe the method and data used in this study in Sect. 2. The mean climatology of the age of air is discussed and compared with observations in Sect. 3. The age variability, the impact of annual cycle and QBO, and the age trend are discussed in Sect. 4. Section 5 provides further discussions and conclusions.

Backward trajectories
Backward deterministic trajectories are calculated using the Lagrangian model TRACZILLA (Legras et al., 2005) which is a modified version of FLEXPART (Stohl et al., 2005). TRACZILLA uses analysed winds to move particles in the horizontal direction and performs direct interpolations from data on hybrid levels. In the vertical direction, it uses either pressure coordinate and Lagrangian pressure tendencies, or potential temperature coordinate and heating rates. In the first case, we denote the trajectories as kinematic and in the second case as diabatic following a convention established by Eluszkiewicz et al. (2000). At each vertical level, particles are initialised over a longitude-latitude grid with 2 • resolution in latitude and an almost uniform spacing in longitude of 2 • / cos(φ), where φ is the latitude, generating 10 255 particles on each level. For convenience the vertical levels of the initial grid are chosen to be the hybrid levels of the ECMWF model. In order to encompass the whole stratosphere at any latitude, the 30 levels from about 400 hPa (varying according to the surface pressure) to 2 hPa are selected. Trajectories starting below the tropopause are immediately stopped and therefore do not induce any computational cost. Particles

Data
The wind data and heating rates used in this study have been produced by the ERA-Interim reanalysis of ECMWF (Dee et al., 2011). This reanalysis uses a 12 h 4D-Var assimilation cycle with a T255 partially desaliased horizontal truncature in spherical harmonics and 60 hybrid levels in the vertical from the surface to 0. Several studies have shown that winds from analysis or reanalysis are noisy and induce unrealistic diffusive transport and too fast apparent Brewer-Dobson circulation in the stratosphere (Schoeberl et al., 2003;Meijer et al., 2004;Scheele et al., 2005). This effect is mostly noticed in the vertical direction where velocities are naturally very small. There are two main reasons for this behaviour. The first is the gravity wave noise induced by the assimilation system. Such noise is transient and dampened during subsequent evolution so that medium-range forecasts exhibit less diffusion than the analysis (Stohl et al., 2005;Legras et al., 2005). This effect is pronounced in assimilation systems using 3D-Var assimilation, like in the ERA-40, and is significantly reduced with 4D-Var assimilation, like in the ERA-Interim. The second reason lies in the fact that archived analysis used for off-line transport studies are instantaneous winds typically sampled at 6 h interval. As a result, fast perturbations with time-scale smaller than 6 h are under-sampled. In the limit of very fast uncorrelated perturbations, a sampling with interval τ induces a spurious diffusion which is proportional to τ . Other reasons might be found in the parameterisations of gravity-wave drag, the representation of convection and the radiative calculations.
The undersampling effect can be reduced by using higher sampling rates at 3 h resolution (Stohl et al., 2005;Legras et al., 2005) or averaging the wind field (Schoeberl et al., 2003;Schoeberl and Dessler, 2011) as both tend to reduce the noise. There are indications (Pisso and Legras, 2008) that increasing the sampling rate to 1 h does not improve the noise in the stratosphere with current generation of reanalysis.
There are several reasons for which the vertical motion, as represented by heating rates in isentropic coordinates, is expected to be less noisy than that represented by vertical velocities in pressure coordinates. Using isentropic coordinates in the vertical separates the fast isentropic motion from slower vertical cross-isentropic motion in the stratosphere and avoids spurious numerical transport effects when particles move with respect to oscillating isobaric surfaces. Another reason is that the heating rates are usually archived as accumulations over finite periods and not as instantaneous values like the velocities. Consequently, the heating rates integrate the noisy fluctuations and are much smoother in time and space than the kinematic velocities.
Data from the ERA-Interim have been used from 1979 to 2010. Since backward trajectory calculations are performed over a duration of 10 yr, the age of air has been estimated monthly over a 22-yr period between 1989 and 2010. For the sake of comparison with the ERA-40 reanalysis and the calculations of Monge-Sanz et al. (2007), some integrations with a perpetual 2000 year have also been performed.

Metric of Brewer-Dobson circulation
The residence time of a particle in the stratosphere since it has crossed the tropopause, is defined as the age of air (Waugh and Hall, 2002) and is a common metric of the Brewer-Dobson circulation. As each air parcel results from the mixing of a large number of particles with different trajectories within the stratosphere, the age is actually distributed over a range of values for all the particles contributing to a given parcel. This distribution is denoted as the age spectrum which can be mathematically defined (Kida, 1983;Hall and Plumb, 1994;Waugh and Hall, 2002) as generated by a Green function describing the probability that a particle located at the tropopause at time t −τ is found within the considered parcel at time t. The first moment of the distribution of τ is the mean age.
The age of air can be retrieved from trace chemical species, such as sulfur hexafluoride SF 6 and carbon dioxide CO 2 which are well-mixed in the troposphere with a known trend and are nearly passive tracers in the stratosphere (SF 6 is only oxidized in the mesosphere and CO 2 has limited sources through the oxydation of CH 4 ) (Andrews et al., 2001a;Waugh and Hall, 2002;Stiller et al., 2008;Garcia et al., 2011). Such quantities have been measured from aircraft and balloons for several decades (Andrews et al., 2001a) and more recently from satellites (Stiller et al., 2008(Stiller et al., , 2012Foucher et al., 2011). These data are used in this study in order to compare modelized ages to observations. The comparison to more recent observations of CO 2 stratospheric profiles from ACE-FTS (Foucher et al., 2009(Foucher et al., , 2011   different life times is available (Andrews et al., 1999(Andrews et al., , 2001bSchoeberl et al., 2005;Bonisch et al., 2009), but in practice only the mean age can be retrieved without ad hoc assumptions on the age distribution.
In our backward calculations, the age along a given backward trajectory is obtained as the time of first crossing of the tropopause, defined as the lower envelope of the surfaces θ = 380 K and |P | = 2 × 10 −6 K kg −1 m 2 s −1 where P is the Ertel potential vorticity. The mean age for a given box in latitude and altitude (typically 2 • × model level spacing) and for a given month is calculated as the average in longitude over all particles falling within this box. Owing to the quasiuniform spread of the discrete trajectories at the initialisation stage, the average is made over 180 particles at the equator and over 67 particles at 68 • N or S. Latitudes closer to the pole are grouped into enlarged latitude bins (69 • -73 • , 73 • -77 • , 77 • -81 • , 81 • -90 • ) to avoid large fluctuations due to the reduced number of particles. Further averaging over time is performed to improve statistics and to reduce noise. These averaging procedures are a simple way to account for mixing in the stratosphere and gather within each box a distribution of particles with different histories.
As observed by Scheele et al. (2005), the number of backward trajectories launched at a given date and remaining within the stratosphere after some delay τ decreases exponentially with τ . Figure 1 shows that this law is indeed very well satisfied for τ > 3 yr with an exponential decrement b = 0.2038 yr −1 for the mean decay and that the standard deviation from the mean (when each month is considered separately) decays at the same rate. After 10 yr of backward motion, 88 % of the particles launched within the stratosphere have met the tropopause. We follow Scheele et al. (2005) in using this property to correct the estimated ages for the truncature of trajectory lengths at 10 yr. If we define F (τ ) as the probability density of the age τ , the mean age is (1) The truncated version of this integral, up to t f = 10 yr, can be calculated explicitly from the trajectory calculations. Assuming that F (τ ) = F (t f ) exp(−b(τ − t f )) for t > t f , the mean age can be estimated as In practice, the calculation is discretized in the following way. For a total of N particles, those with ages under t f , which have crossed the tropopause during backward integration before t f , are distributed in K ages bins between 0 and t f with n i particles in the bin i centered on time t i . If M f particles remain within the stratosphere at time t f , the corrected mean over this ensemble is The decrement coefficient b, shown in Fig. 1 as a time average, has been calculated over the whole stratosphere for each month. It varies by ±0.02 yr −1 over time. It has also been calculated as a 22-yr mean for each latitude and altitude box. The resulting correction to the mean age varies from zero to almost two years at high altitude and latitude. The impact of choosing one definition of the decrement or the other does not change the estimated age by more than 3 %. Hence, the correction which is not negligible per se is quite insensitive to the arbitrary details of the calculations. Notice, however, that the mean value of b differs from our value by almost a factor 2 in Li et al. (2012) who found b ≈ 0.36 using the GEOSCCM model. The comparison with this calculation is further discussed below.

Global distribution of the mean age
The mean diabatic age of air, obtained with diabatic trajectories, is calculated as a function of latitude and altitude after averaging over the 22-yr dataset between 1989 and 2010. The left panel of Fig. 2 shows that mean age contours follow the tropopause except in the tropics where they bend up as a result of the tropical upwelling. Gradients of the age of air are concentrated within the extra-tropical lowermost stratosphere with approximately 0.5 yr per km. In the mid extratropical stratosphere, the mean age of air varies between 6 and 7.5 yr, with maximum values near the poles, and is older in the Northern Hemisphere above 25 km than in the Southern Hemisphere. The tropical pipe (Neu and Plumb, 1999), which confines the ascending branch of the Brewer-Dobson circulation above about 22 km, is revealed by young air moving upward in the tropics. This tropical pipe is slightly shifted from the equator with a maximum near 5 • S. Its relative isolation is visualised by the horizontal age gradients on its northern and southern edges. The right panel of Fig. 2 shows the standard deviation of the mean diabatic age calculated using the equivalent sampling size (see Appendix) that accounts for the time correlation of monthly ages. This standard deviation shows that the mean diabatic age is estimated with good accuracy within the framework of the ERA-Interim, with patterns that clearly offset the level of fluctuations. The variability of the mean age concentrates within a limited range of altitudes between 20 and 30 km and, as we shall see below, that the maximum of variance is correlated with the maximum of QBO wind modulation.
The difference between kinematic and diabatic ages, see Fig. 3, has been calculated over 35 months within the period 2006-2009. The pattern is quite unexpected with kinematic ages being older than diabatic ages in the lower stratosphere of the Southern Hemisphere and being younger between 22 and 40 km in the Northern Hemisphere. There is almost no difference in the southern mid-stratosphere and only a thin layer of older kinematic ages is visible in the northern lowermost stratosphere. It is usually found that kinematic velocities are noisier than diabatic heating rates resulting in a spurious vertical diffusion and a bias towards younger air (Schoeberl et al., 2003;Ploeger et al., 2010). However, the pattern of Fig. 3 cannot be explained by such a simple argument. In order to understand better the relation between air parcel origins and ages, Fig. 4 shows the distributions of maximum vertical excursion and of altitude of tropopause crossing. The distributions are shown separately for particles launched below 113 hPa in the extra-tropical lowermost stratosphere (lower row) and those launched at this level and above in the tropics and the extra-tropics (upper row), in the region of the stratosphere denoted as the overworld (Holton et al., 1995). For overworld parcels, Fig. 4 shows that most of the entries to the stratosphere occur through the tropical tropopause between the isentropic levels 370-380 K. The histogram of maximum vertical excursion shows three maxima. The first one below 500 K corresponds to the fast branch of the Brewer-Dobson circulation which is bound to the lower stratosphere. The two other maxima are associated with the tropical pipe. The plume of air rising through the pipe is progressively stripped by detrainment to the mid-latitudes. Most particles reach a maximum value under 1500 K. Above this level, there is very little leakage from the tropical pipe between 1800 K and 2500 K and the third maximum near 2800 K is associated with particles reaching the top of the mesosphere in the model.
For particles initialised in the extra-tropical lowermost stratosphere, Fig. 4 shows that the majority of tropopause crossings still occur between 370 and 380 K in the tropics but a significant proportion of particles enter the stratosphere through the subtropical and extra-tropical tropopause at lower potential temperatures down to 300 K. The maximum vertical excursion is mainly contained within the 300-450 K range with a peak at 380 K. Only a small portion of the particles (about 6%, not shown), have maximum vertical excursion exceeding 700 K.
Hence, the low value and the strong gradient of the mean age above the extra-tropical tropopause are due to the combined effect of isentropic mixing of tropical and extra-tropical air across the subtropical tropopause and the fast shallow branch of the Brewer-Dobson circulation (Hoor et al., 2004;Bonisch et al., 2009). Although the deep branch of the Brewer-Dobson is important for the distribution of ages in the stratosphere and for stratospheric chemistry, it processes only a small portion of the air which circulates within the stratosphere and the air found within the lowermost extra-tropical stratosphere (except within the winter polar vortex) has been mainly processed through the shallow branch.

Comparison with observations and models
As a basis for comparison, we use the age of air obtained from in situ aircraft data prior to 1998 reported in Andrews et al. (2001a) and the ages derived from MIPAS retrieval of SF 6 in 2002(Stiller et al., 2008. Figure 5 shows that these two estimates overlap in the mid-latitudes but the SF 6 ages are older at high latitudes. This is consistent with the impact of photochemical dissociation of SF 6 in the mesosphere which contaminates the stratospheric air within the winter polar vortex (Waugh and Hall, 2002). However, according Annual Mean Age in ERA−Interim for Kinematic (Z) and Diabatic (Dia) trajectories runs at 54.62hPa SF6 MIPAS: Stiller & al., 2009 Variance of Mean age−of−air In−situ CO2: Andrews & al., 2001 In−situ SF6: Elkins & al., 1996 Annual  (Elkins et al., 1996;Waugh and Hall, 2002). (Dots and error bars): compilation of mean ages from airborne observations of CO 2 until 1998 by Andrews et al. (2001a). The error bar shows the statistical uncertainty of the mean age. The statistical uncertainty of the mean diabatic age is shown as ± one standard deviation (cyan shaded area) but this uncertainty is small enough to be hidden by the thickness of the red curve.
to Stiller et al. (2008), the systematic errors of age retrieval from SF 6 are such that ages are 0 to 0.5 yr too young in the lower stratosphere, even if they remain older than other observations. Figure 5 shows that the ERA-Interim mean diabatic ages for the period 1989-2010 (red curve) are in good agreement with the aircraft observations except at high latitude. They are generally smaller than the MIPAS SF 6 ages by about 1 yr except in the southern mid-latitudes where the agreement between observations and simulation is the best. Consistently, the observations are less dispersed in this region. In the tropics the SF 6 MIPAS mean ages are about twice that of the ERA-Interim and in situ observations of SF 6 and CO 2 . These comparisons should be appreciated with the reservation that observed and simulated ages are obtained over overlapping, albeit non identical, periods.
As already noticed the kinematic trajectories tend to produce significantly older ages in the Southern Hemisphere (black curve) at 20 km. These kinematic trajectories have been calculated for two years only, 2007 and 2008, but the discrepancy is meaningful because the diabatic ages averaged over the same years (blue curve) do not depart significantly from the 22-yr mean. This is contrasted with the kinematic trajectories calculated with winds from the ERA-40 reanalysis which tend to systematically produce younger ages (see Monge-Sanz et al., 2007, 2012, and cyan curve in Fig. 5). In the ERA-Interim, the statistics produced by kinematic and diabatic trajectories differ less than in the ERA-40 (Liu et al., 2010) but the difference is reversed. Schoeberl and Dessler (2011), using the MERRA reanalysis, found kinematic ages older than diabatic ages in the whole stratosphere and at the same time found an excessive vertical diffusion associated with the kinematic trajectories. Why kinematic trajectories produce longer residence time in the stratosphere or the type of pattern seen in In situ SF6: Ray & al., 1999 In situ CO2: Andrews & al., 2001& Boering & al., 1996 Air samples ouside vortex SF6: Harnisch & al., 1996 Air samples inside vortex SF6: Harnisch & al., 1996SF6 MIPAS: Annual Mean Profile Stiller & al., 2011 Vertical profiles of mean age of air in ERA−Interim from diabatic heating rate  (Boering et al., 1996;Andrews et al., 2001a), SF 6 (triangle) (Ray et al., 1999) and whole air samples of SF 6 (square outside vortex and asterisk inside vortex) (Harnisch et al., 1996). produce younger ages, and reinforced exchanges between the tropics and the mid-latitudes should produce older ages in the tropics and younger ages in the mid-latitudes. None of these patterns is observed in Schoeberl and Dessler (2011) or in our study.
We stress that our calculations are all based on full historical records of velocity fields and heating rates over the length of the integration. In a number of previous studies, simulations using perpetual repetition of a given year have been used. This choice leads to considerable fluctuations of the age of air. For instance, the kinematic ages obtained for ERA-Interim based on a 2000 perpetual are significantly older than the 22-yr average (see Fig. 5) while its ERA-40 counterpart provides much too young ages. Large fluctuations, positive or negative, are also observed for diabatic trajectories calculated over perpetual years for 2000 and other years (not shown).
When compared with Chemistry-Climate-Model (CCM) estimates , our ages based on Lagrangian trajectories are usually older, by about 1 yr in the tropics above 30 km and often 2 yr at mid and high latitudes. The horizontal gradient between the tropical pipe and the mid-latitudes is also stronger. However, a detailed comparison with a recent study of the age of air in the GEOSCCM  reveals that the patterns of the age of air distributions, including its annual variations, are strikingly similar, even if the shift in the mean age just mentioned is still oberved. The comparison is performed in the Supplement where our results have been redrawn to produce figures that can be directly compared with those of Li et al. (2012). The differences in the mean ages may be partially due to differences in the numerical representation of tracer advection (Eluszkiewicz et al., 2000) as non diffusive Lagrangian calculations tend to produce older ages than other methods. Quite interestingly, it is shown in the supplement that discarding all particles travelling above θ = 1800 K does not change the patterns of ages distribution but improves considerably the agreement with GEOSCCM ages. This suggests that trapping of particles near the lid of the model in the ERA-Interim might lead to an old bias.

Vertical profiles of the mean age
A detailed comparison of the vertical profiles of the calculated mean ages with those derived from observations of middle stratosphere balloon flights (Andrews et al., 2001b;Ray et al., 1999) and from SF 6 MIPAS profiles is shown in Fig. 6. According to Stiller et al. (2008), the systematic errors of SF 6 ages are such that ages are 0 to 1 yr too old between 25 and 35 km. In the tropics, the ages from SF 6 MIPAS are higher than those from in situ measurements at all altitudes. The diabatic ages are in good agreement with the in situ measurements up to 28 km and with SF6 ages above 30 km. The Atmos. Chem. Phys., 12, 12133-12154, 2012 www.atmos-chem-phys.net/12/12133/2012/ seasonal dispersion is much smaller than the discrepancy between in situ and satellite data. The mean diabatic age increases almost uniformly in z from 18 to 34 km at a rate of 0.35 yr km −1 .
In the mid-latitudes, there is a good agreement between in situ and satellite observations, except for the SF 6 profiles of Harnisch et al. (1996) above 25 km. The diabatic ages follow the main group of observations, being slightly on the younger side during summer below 20km. The age increases from 16 to 28 km at a rate of 0.4 yr km −1 and exhibits only a weak vertical gradient above 28 km, consistently with the findings by Waugh and Hall (2002).
At high latitude during winter, the descent of mesospheric air within the vortex and on its edge induces strong contrast between in vortex and out of vortex air which shows up as the difference between winter and summer profiles below 22 km. Above 25 km there is a small difference between winter and summer as in the mid-latitudes.
The mean diabatic ages agree with SF 6 ages derived from Harnisch et al. (1996) in vortex observations and MIPAS up to 28 km. They depart from other in situ observations by one to two years over the whole altitude range.
Above 28 km, the SF 6 ages from MIPAS exhibit a positive departure which is consistent with the partial photo-chemical dissociation of SF 6 in the mesosphere (Stiller et al., 2012).
The  -Sanz et al. (2012) are two folds: a) they use a CTM when we are using Lagrangian calculations and b) they use 2000 perpetuals when we use full historical records for all calculations. As already mentioned, pure Lagrangian calculations are not introducing artificial numerical diffusion which tends to limit old ages. We have some reservation about using single year perpetuals as there is a risk of introducing a bias by freezing the QBO oscillation in one of its phases.
The whole comparison suggests, however, that the ERA-Interim tend to produce older ages than observed, especially at high latitudes.

Age spectrum
As a significant portion of particles remain in the stratosphere with old ages (see Fig. 1), it is important to consider not only the mean age but also the age spectrum. Figure 7 shows that there is a clear distinction between the tropical and the extratropical spectra at 20 km. In the tropics, the distribution of ages is mono-disperse and compact, and decays rapidly to zero for ages above 1 yr indicating that very few particles with old ages return to the tropics from mid-latitudes. In the extra-tropics, the peak is at about 0.5 yr, which is small compared to the mean corrected age of about 3.5 yr. The ages exhibit a long flat tail which extends well to large ages. This age distribution of ages corroborates the existence of fast and slow branches of the Brewer circulation (Bonisch et al., 2009(Bonisch et al., , 2011 even if a secondary maximum is not seen. The fast branch is associated with the particles which have travelled directly and rapidly from the tropics to the mid-latitudes through quasi-isentropic motion (Haynes and Shuckburgh, 2000a;Hoor et al., 2005;Shuckburgh et al., 2009), staying at levels below 450 K. The slow branch corresponds to the deep Brewer-Dobson circulation in which the particles enter the tropical pipe and circulate to high altitudes in the stratosphere. A strong seasonal modulation of the fast branch is observed in the Northern Hemisphere with a younger peak during summer than during winter. The modulation is smaller in the Southern Hemisphere. This variation is associated with the seasonal modulation of the subtropical jet and the meridional exchanges which is larger in the Northern Hemisphere than in the Southern Hemisphere (Hoor et al., 2005;Shuckburgh et al., 2009).
At higher altitudes (see Fig. 7), the age distribution in both the tropics and the extra-tropics shifts to older ages. The tail of the tropical distribution gets thicker with increasing altitude but remains much less developed than the extra-tropical tail up to 30 km, in agreement with the relative isolation of the tropical pipe (Neu and Plumb, 1999). The peak of the tropical distribution is 1 to 2 yr younger than the broader main maximum of the extra-tropical distribution. A striking feature is the presence of oscillations in the distributions with an interval of one year between maxima. These oscillations, already mentioned by Reithmeier et al. (2008) and Li et al. (2012), reach a maximum amplitude at the modal age. They propagate towards old ages with the seasonal cycle and it is remarkable that the phase of the oscillation is the same at all altitudes and in both hemispheres. During boreal winters the maximum occurs for integer ages; a positive shift of three months is added for spring and so on for the other seasons. This is consistent with the interpretation of Reithmeier et al. (2008) that the oscillations are entirely due to the modulation of the mass flux entering the stratosphere which is indeed maximum during boreal winter (Seviour et al., 2011). The decay of amplitude with the age, which is almost perfectly exponential in Li et al. (2012) is due to the repeated and multiplicative action of mixing. A detailed comparison of our results with those of Li et al. (2012) is provided in the supplement http://www.atmos-chem-phys.net/12/12133/ 2012/acp-12-12133-2012-supplement.pdf. We find again an excellent agreement in the patterns of the age spectrum and its seasonal variability with the sole exception that our phase modulation is still in phase with a maximum winter upwelling at high latitude, in agreement with  Li et al. (2012) find that air leaving the tropical region in summer has more chance to reach the polar stratosphere. Finally, it should be stressed that the seasonal modulation of the spectrum does not necessarily show up as a modulation of the mean age when the age distribution is fairly flat as it occurs at 31 and 35 km. This can be checked in Fig. 6 where the profiles are very close in winter, summer and in the annual mean at such levels.

Variability and trends
The age of air contains an integrated footprint of the variability of the Brewer-Dobson circulation. The first hint on variability is provided by considering the temporal variation of the number of remaining particles with respect to the mean decay shown in Fig. 1. It is visible (upper left panel of Fig. 8) that the variations are dominated by the annual cycle. The mean annual cycle (lower left panel) exhibits a negative deviation during winter and a positive deviation during summer, implying that more particles cross the tropopause during boreal winter than during summer, consistently with the winter intensification of the tropical upwelling. This annual cycle propagates through ages and it is useful to notice that the amplitude of the cycle is largest for ages of 2 to 3 yr which indicates the duration over which the exiting particles strongly feel the annual cycle. The amplitude decays as age gets older since and remaining particles get uniformly distributed within the stratosphere. This is consistent with the decaying oscillations observed in the age spectra.
The right panel of Fig. 8 shows the percentage of remaining particles after removal of the mean and of the annual cycle. Special events affecting transfers at the tropopause are seen as discontinuities in the vertical whereas variations of the Brewer-Dobson circulation are seen as the oblique patterns.
The most prominent feature is clearly associated with the Pinatubo eruption in June 1991 which induces a reduction of the tropopause crossings and a slowing of the Brewer-Dobson circulation. The impact is strongest on particles with ages of 2 to 3 yr but extends to much older ages and over most of the following decade.
The second smaller pattern after 2008 is possibly due to the cumulative effects of the eruptions of the Soufrières Hills Atmos. Chem. Phys., 12, 12133-12154

Regression method
The temporal evolution of monthly mean diabatic ages at specific altitudes and latitudes (bined as described in Sect. 2.3) has been analysed using a linear response model over the 22 yr of data available from TRACZILLA integrations. This model yields where qbo is a normalised quasi-biennial oscillation index from CDAS/Reanalysis zonally averaged winds at 30 hPa and enso is the normalised Multivariate El Niño Southern Oscillation Index (MEI) Timlin, 1993, 1998), both provided by the NOAA website. The coefficients are a linear trend a, the annual cycle C(t) (12 coefficients), the amplitude b 1 and the delay τ qbo associated to QBO and the amplitude b 2 and the delay τ enso associated to ENSO. The constraint applied to determine the 17 parameters a, b 1 , b 2 , τ qbo , τ enso and C is to minimise the residual ǫ(t) in the least square sense. As the combination of amplitude and delay introduces a non linear dependency, there are multiple minima which are sorted out as described in Sect. 4.3.
Here we discard the influence of Pinatubo because there is no simple index suitable to describe the effect of this single event on the Brewer-Dobson circulation and we also neglect solar forcing, because our dataset covers only two solar periods.

Annual cycle
The annual cycles calculated by the minimising procedure (which takes in account all the factors of variability together) or by a simple monthly composite over the 22 yr turn out to be almost identical. Figure 9 shows the amplitude of the annual cycle and its phase calculated by fitting a pure annual cosine variation to the full annual cycle. The phase of the cosine is defined to be zero in mid-January. The amplitude is maximum in the extratropical lowermost stratosphere with peak values of about one year at all latitudes in the Northern Hemisphere and of half a year, except at high latitudes, in the Southern Hemisphere. The phase is in opposition between the two hemispheres, the maximum being in mid-March for the Northern Hemisphere and mid-September in the Southern Hemisphere. This maximum signal at the end of the winter is consistent with the reinforced barrier effect of the jet and a stronger descent of old air due to the intensification of the deep Brewer-Dobson circulation during winter (Holton et al., 1995;Waugh and Hall, 2002). In turn, younger ages are observed during summer and autumn. Above 25 km, where the amplitude modulation is less than 0.5 yr, the phase is in opposition with the extra-tropical lowermost stratosphere. At these altitudes, the stronger descent during winter favours the replacement of old air by younger air detrained from the tropical pipe. The maximum modulation in the Southern Hemisphere is found at 60 • S and is associated with the strong descent occurring on the polar vortex edge during late winter and spring (Mariotti et al., 2000).
The enhanced penetration of subtropical and tropical air masses in the extra-tropics during summer and autumn, due to a decreased barrier against quasi-horizontal transport Atmos. Chem. Phys., 12, 12133-12154, 2012 www.atmos-chem-phys.net/12/12133/2012/ and mixing at the subtropical jet, is consistent with previous works on stratosphere-troposphere exchanges, based on tracer's observations (Hoor et al., 2005;Krebsbach et al., 2006;Sawa et al., 2008) or on model simulations (Chen, 1995;Haynes and Shuckburgh, 2000b;Sprenger and Wernli, 2003). It is also consistent with our findings on the age spectrum and those by Bonisch et al. (2009). In the tropics, the annual modulation above the tropopause and in the tropical pipe rarely exceeds 0.5 yr. In the Northern Hemisphere, it is confined below 30 km with a maximum during summer. In the Southern Hemisphere, it is confined between 30 and 40 km with a maximum at the end of the winter. As already mentioned, the annual modulation is very similar to that described by Li et al. (2012) The seasonal modulation is shown with more details at 55 hPa (about 20 km) in Fig. 10 where the diabatic ages are compared with ages derived from MIPAS SF 6 data (Stiller et al., 2012). It is visible that the Northern Hemisphere modulation is larger in the tropics and in the extra-tropics while the Southern Hemisphere modulation is larger in the subtropics. At high latitudes, both hemispheres exhibit the same amplitude. It is also visible that if SF 6 ages and diabatic ages are of similar amplitude at mid and high latitudes, their seasonal variation are not related.

Quasi-Biennal-Oscillation and ENSO
Because of the presence of lags in the QBO and ENSO terms in Eq. (4), the problem is non linear and the residual may have multiple minima as a function of the parameters. In order to determine the optimal values of τ qbo and τ enso , the residual is first minimized at fixed lag and then over a range of lags. This is done in sequence for QBO and ENSO. Figure 11 shows the variations of the QBO amplitude coefficient and the residual amplitude as a function of latitude and lag at several levels in the vertical, roughly corresponding to 20, 25 and 30 km height. In most cases, the minimum residual corresponds also to a maximum of the QBO amplitude coefficient in absolute value and the QBO correlates with the age over more than one period. The optimal lag strongly depends on latitude, varying, e.g., by more than a year between 0 and 20 • S at 25 km. Figure 12 shows the amplitudes and lags of the QBO and ENSO contributions to the variability of the age of air. The QBO modulation reaches 0.6 yr with a lag of about 8 months in the tropics near 30 km, with a stronger component in the Northern Hemisphere. This region has already been identified as displaying the largest variability in the age of air in Fig. 2. The influence of the QBO extends upward in the tropical pipe and towards the extratropical stratosphere in both hemispheres with amplitudes of the order of 3 months. The phase lag with respect to the wind at 30 hPa is fairly symmetric with respect of the equator and varies most rapidly with latitude near 25 km. The dependence on ENSO is much less pronounced (less than 0.2 yr) and bound mostly to the lower stratosphere in the Northern Hemisphere. The standard deviation of the residual in Eq. (4), shown in the right panel of Fig. 12, is larger than the signal explained by QBO and ENSO and is maximum in the same regions as the QBO. Hence, the variability not linked to QBO and ENSO dominates the age of air at any location in the stratosphere.

Trends
One should be cautious when estimating a trend in the ERA-Interim, because the reanalysis system cannot be considered time invariant, due to numerous changes in the observations, in particular the introduction of new satellite instruments. Such changes are liable to induce biases in the atmospheric circulation in spite of the attention devoted to avoid them (Dee et al., 2011). It is nevertheless useful to determine the trends, and to compare and eventually reconcile them with observations.
As for the annual cycle, the trends calculated by the minimising procedure (which takes in account all the factors of variability together) or by a simple linear fit turn out to be almost identical. The trend is shown as a function of latitude and altitude in Fig. 13 (left panel). The trend is negative within the lower stratosphere with a larger magnitude in the tropics and the Southern Hemisphere, of the order of −0.3 yr dec −1 , than in the Northern Hemisphere. The trend is positive in the extra-tropics above 25 km. In the tropics, the situation is contrasted between the Southern Hemisphere where the trend is negative up to 33 km and the Northern Hemisphere where it is positive above 28 km. The maxima of the trend are located where the amplitude of the QBO modulation is also the largest.
The significance of the trend has been assessed, following von Storch and Zwiers (1999), by performing a Student's ttest among the 264 months of our record using an equivalent number of degrees of freedom calculated as in Zwiers and von Storch (1995) and Bence (1995) (see Appendix). This equivalent number ranges between 35 and 80 in the region of maximum negative trend. The right panel of Fig. 13 shows the one-sided p-value of the Student's test for the hypothesis of a null trend. It is visible that the whole region of large negative trend is highly significant. Although such a simple test is known to overestimate the significance in many cases (Zwiers and von Storch, 1995) and we neglect sub-monthly contribution to the mean age variance, this is an indication that the negative trend in the lower stratosphere is a robust Atmos. Chem. Phys., 12, 12133-12154 feature. On the opposite, we do not find any area of consistent significance above 25 km except perhaps near 60 • N and 30 km. In the tropics near 30 km, the p-value is 0.30. The observed negative trend in the lower stratosphere is consistent with the finding of an acceleration of the shallow branch of the Brewer-Dobson circulation by Bonisch et al. (2011). The positive trend in the mid and high latitudes above 25 km, of the order of 0.2 yr dec −1 is consistent with the findings of Engel et al. (2009) andStiller et al. (2012). It is, however ten times larger than the value found by Monge-Sanz et al. (2012) using the same data. In order to assess better the contributions to the age of air, Fig. 14 shows three illustrative cases. The first case is chosen at 20 km and 60 • N, because the amplitude of the annual cycle is the largest and indeed dominates the variability of the age of air, the QBO component being much smaller. The second case is chosen at 28 km and 4 • N, because the QBO component dominates the variability over the annual cycle. The ENSO component makes a negligible contribution in both cases. The trend is negative of −0.29 yr dec −1 in the first case and slightly positive of 0.11 yr dec −1 in the second case. The age fluctuations are large in the second case, so that the trend is hardly significant. In the first case, it is visible that the last section of the record, after 2004, would suggest a positive trend while the negative trend over the whole period has a p-value of 0.02. The third case is chosen at 19.7 km and 40 • S, because the trend exhibits its largest value of −0.56 yr dec −1 . This is obtained in spite of a flattening or a very slight increase over the last period but it is also visible that the period following the Pinatubo eruption exhibits increased ages which significantly contribute to the trend. The amplitude of this effect is tested by recalculating the trend after removing the period mid-1991 to the end of 1994, leading to a reduced value of −0.33 yr dec −1 . It is thus clear that the age of air undergoes decadal scale variability, induced by volcanoes and other less known processes, which may affect trends calculated over durations of a few decades. Our own estimate, based on 22 yr of estimated ages is prone to such effect and our statistical test is not made against such decadal variability.
We compare now our calculations to the age of air obtained by Stiller et al. (2012) using SF 6 data from the MIPAS instrument on board ENVISAT over the period 2002-2010 (see Fig. 15). The comparisons are performed at three altitudes (16.5, 20 and 25 km) and four bands of latitudes (30 • S-20 • S, 20 • S-10 • S, 40 • N-50 • N, 70 • N-80 • N). There is an old bias of SF 6 ages at 16.5 km at all latitudes and at 20 km in the tropics with respect to our calculations. The tropical bias, which was already noticed on Fig. 5 appears as "robust" in the sense that the dispersion of SF 6 ages is much smaller than the bias. The other panels show a good agreement between our estimations and SF 6 ages. A striking feature is that the age of air tends to flatten or to rise in most locations over the 2002-2010 period, even when a clear negative trend is found over the 22-yr period. This is also in agreement with Stiller et al. (2012) who found a well-spread positive trend over 2002-2010 although there is a disagreement in some locations, for example at 25 km and 70 • N-80 • N. This is, however, where the statistics are noisy and the trends unreliable.

Conclusions
We have performed direct calculations of the stratospheric age of air using Lagrangian trajectories guided by reanalysed winds and heating rates from the ERA-Interim. This study is based on 32 yr of data and provides estimates of the age over the last 22 yr of the dataset . This study complements previous works on the Brewer-Dobson circulation in the ERA-Interim (Iwasaki et al., 2009;Garny et al., 2011;   Our analysis corroborates the significant improvement of the stratospheric circulation in the ERA-Interim compared to the older ERA-40, as previously noted by other recent works (Monge-Sanz et al., 2007, 2012Fueglistaler et al., 2009b;Dee et al., 2011;Seviour et al., 2011). In contrast with the ERA-40, diabatic heating rates provide younger ages than the kinematic velocity, a feature shared with the MERRA reanalysis (Schoeberl and Dessler, 2011) but which is not yet properly understood.
On the overall, the agreement between the diabatic ages and observations is very good at low and mid-latitudes but for the SF 6 ages in the tropics. At high latitude, the diabatic ages agree with SF 6 ages but are older than other observations. The comparison with GEOSCCM ) (see Supplement) demonstrates that a state-of-the-art CCM and the ERA-Interim generate very similar patterns of the stratospheric age, in the mean, the annual variations and the age spectrum. The old bias of the ERA-Interim ages with respect to GEOSCCM and other CCMs  entirely disappears, without changing the patterns, when parcels reaching above θ = 1800 K are discarded from the calculation. This is puzzling because CCMs have very different upper lids, sometimes above the ERA-Interim (like GEOSCCM) but many times below, and suggests that the old bias of ERA-Interim with respect to the CCMs may arise from a variety of factors.
The age spectrum corroborates that the lower extratropical stratosphere contains young air coming from the mass residual circulation (Birner and Bonisch, 2011) and horizontal mixing, both combined in the lower branch of the Brewer-Dobson circulation. At higher altitude, the age spectrum reveals peaks spaced by annual intervals that propagate on the spectrum with the seasonal cycle and are associated with the annual modulation of the Brewer-Dobson circulation. These modulations do not necessarily show up or weakly when the sole mean age is considered. While the age spectrum is hardly accessible from observations without hypothesis on its shape (Andrews et al., 2001b;Schoeberl et al., 2005;Bonisch et al., 2009), we suggest that it would be useful to expand its usage within the framework of model inter-comparisons as it contains much more informations on the stratospheric circulation than mean age only.
Atmos. Chem. Phys., 12, 12133-12154  The variability of mean age as a function of latitude and altitude had been analysed as a linear combination of several contributions: annual cycle, QBO, ENSO and trend. The annual variability dominates in the extra-tropical lower stratosphere below 25 km with higher amplitude, reaching 1.5 yr, in the extra-tropical Northern Hemisphere than in the Southern Hemisphere. The age reaches its maximum during February-March in the Northern Hemisphere and during August-September in the Southern Hemisphere. This is in agreement with a reduction of the shallow branch of the Brewer-Dobson circulation and an intensification of the deep branch during the winter in each hemisphere which has been described in many studies (e.g. Iwasaki et al., 2009;Garcia et al., 2011;Seviour et al., 2011) but it has little relation with the annual modulation as derived by Stiller et al. (2012). In the tropics, the amplitude of the annual modulation does not exceed six months and is basically in opposition with the modulation of the tropical upwelling which reaches its maximum during boreal winter.
The QBO modulation is mostly pronounced within the tropics between 25 and 35 km where its maximum reaches 6 months. The ENSO signal is found to be small and noisy, with a maximum value of about 3 months, and limited to the lower northern stratosphere. This does not contradict the recent finding that warm ENSO events accelerate the upwelling in the lowermost tropical stratosphere (Calvo et al., 2010).
According to these ERA-Interim calculations, a negative trend of the order of −0.2 to −0.4 yr dec −1 is found in the lower stratosphere of the Southern Hemisphere and under 40 • N in the Northern Hemisphere below 25 km. A positive trend of 0.2 to 0.3 yr dec −1 is found in the mid extra-tropical stratosphere above 25 km. The negative trend is significative with respect to a test that ignores the decadal variations of the stratospheric circulation. It is, however, influenced by the Pinatubo eruption and reduced when the Pinatubo years are removed. The positive trend is marginally significant but is consistent with studies based on in situ and satellite observations Stiller et al., 2012). The positive trend in our calculation is increased if only the decade 2000-2010 is considered and this is consistent with the strong positive trend found by Stiller et al. (2012) over the same period.
The whole pattern suggests that the shallow and deep branches of the Brewer-Dobson circulation have not recently evolved in the same direction. The trends are consistent with the recent conclusions of Ray et al. (2010), that the increase in the tropical upwelling documented by Randel et al. (2006) is associated with an acceleration of the shallow branch of the Brewer-Dobson circulation (Bonisch et al., 2011) while the deep branch does not change or slightly weakens. Such pattern is plausible if the shallow and deep branches are controlled by different processes (Gerber, 2012).
CCMs and studies based on residual circulation all predict an increased Brewer-Dobson circulation in the whole stratosphere of 2 to 3 % per decade (McLandress and Shepherd, 2009;Iwasaki et al., 2009;Butchart et al., 2010;Garny et al., 2011;Garcia et al., 2011;Seviour et al., 2011). The discrepancy with our results and those of Engel et al. (2009) andStiller et al. (2012) might arise from the fact that the age of air does not only depend on the residual circulation but also on the stirring and mixing properties of the stratospheric motion. Owing to their limited resolution, CCMs may not necessarily represent well the mixing properties. A separate diagnostic of the isentropic mixing was performed by Ray et al. (2010) who found a large dispersion between reanalysis but concluded that mixing properties may have varied recently. Such study remains to be done for the ERA-Interim.
It is difficult to estimate a reliable long-term trend from only 22 yr of data. There are evidences that volcanic eruptions, ENSO and other badly represented processes like solar variations may affect the stratospheric circulation over decadal scale. An other factor is ozone depletion to which a significant part of the recent evolution of the Brewer-Dobson circulation has been attributed by Li et al. (2008).
There are also limitations due to the fact that a reanalysis dataset like the ERA-Interim is based on a single model but on a constantly changing observation systems. For instance it is known that the introduction of AMSU-A data in 1998 (Dee and Uppala, 2009) or the introduction of radiooccultation data at the end of 2006 (Poli et al., 2010) had significant impact on the stratospheric temperature. To which extend these changing biases may affect the age of air calculations remains to be determined. It is also worth mentioning that the budget of a reanalysis system, both for heat and momentum, is not closed but is biased by the assimilation increment. Such biases may affect trajectory calculations and residual velocities. It is known, however, that the biases are significantly reduced in the ERA-Interim with respect to previous reanalysis (Fueglistaler et al., 2009b;Dee et al., 2011).
Finally the understanding of stratospheric circulation is limited by the poor accuracy and sparseness of available measurements of the mean age of air from tracers. More observations and in situ measurements combined with new processing of available data (see, e.g. Foucher et al., 2011) will improve our understanding.

Appendix A Equivalent sample size
The equivalent sample size for a series of the ages {a i }, where i ∈ [1, N ], is smaller than N because the ages are correlated. We use the Prats-Winstein formula for the equivalent sample size (Bence, 1995;von Storch and Zwiers, 1999): This formula strictly applies to an AR(1) process. By neglecting the periodicities of the age, it provides a conservative lower estimate of the number of degrees of freedom in the series.
The equivalent sample can be used to estimate the standard deviation of the mean as whereā is the mean. It is also used in the Student test for the mean trend.