Retrieving the age of air spectrum from tracers : principle and method

Surface-emitted tracers with different dependencies on transit time (e.g., due to chemical loss or time-dependent boundary conditions) carry independent pieces of information on the age of air spectrum (the distribution of transit times from the surface). This paper investigates how and to what extent knowledge of tracer concentrations can be used to retrieve the age spectrum. Since the tracers considered depend linearly on the transit time distribution, the question posed can be formulated as 5 a linear inverse problem of small dimension. An inversion methodology is introduced, which does not require any assumptions regarding the shape of the spectrum. The performance of the approach is first evaluated on a constructed set of artificial radioactive tracers derived from idealized spectra. Hereafter, the inversion method is applied to model output. The latter experiment highlights the limits of inversions using only parent radioactive tracers: they are unable to retrieve fine scale structures such as the annual cycle. Improvements can be achieved by including daughter decaying tracers and tracers with an annual cycle at the 10 surface. This study demonstrates the feasibility of retrieving the age spectrum from tracers, and has implications for transport diagnosis in models and observations.


Introduction
The transport of surface-emitted tracers strongly influences the composition and chemistry of the atmosphere, as well as the global radiative balance (Riese et al., 2012).In turn, radiatively active species affect the diabatic budget, eventually reshaping the circulation and thus the transport itself.For instance, climate models predict a strengthening of the stratospheric Brewer-Dobson circulation caused by increasing anthropogenic greenhouse gas emissions at the surface (Butchart et al., 2010).
To characterize transport from the surface to a given region of the atmosphere, a number of observational (e.g.Engel et al., 2009) and modeling studies have focused on the average transit time, the mean age of air.However, it has long been acknowledged that the description of transport provided by the mean age is incomplete (e.g.Hall and Plumb, 1994).Large and small-scale turbulent motions lead to mixing, so that a given air parcel is a mixture of air masses with different paths and transit time from the surface (Waugh and Hall, 2002).Strictly, there is a probability distribution of transit time scales for each air parcel, which is known in the stratospheric literature as the age spectrum (Hall and Plumb, 1994;Waugh and Hall, 2002), while the tropospheric literature more frequently uses the acronym TTD for Transit Time Distribution, Holzer et al. (2003).Considering the full age spectrum rather than the mean age allows to separate between different transit times related to different pathways of transport and to disentangle their potentially contrasted evolutions with climate change (see Ploeger and Birner, 2016, and references therein).It also enables to understand the air composition in a number of species without restricting to inert, linearly increasing tracers (Schoeberl et al., 2000).
Up to date, the stratospheric age spectrum has mainly been estimated in models, using either Lagrangian trajectories (e.g Reithmeier et al., 2008;Diallo et al., 2012) or a set of artificial pulse tracers initialized in the lowest model layer (Li et al., 2012;Ploeger and Birner, 2016).Only a handful of studies (Andrews et al., 1999;Johnson et al., 1999;Schoeberl et al., 2005) have attempted to infer the age spectrum from observed tracers, and all assumed either a given shape for the distribution or steadiness of the flow.Many tracers, however, bear the imprint of specific regions of the age spectrum (e.g.Waugh et al., 2003Waugh et al., , 2013;;Orbe et al., 2016) and, combined together, may provide information on the entire transit time distribution.
In this study, we propose a general methodology for retrieving the age spectrum from the concentrations of more general (non-pulse) tracers, which may undergo chemistry and have time-dependent sources.The basic idea is to consider the tracer contents as the images of the age spectrum through a known forward model, and to pose the retrieval of the age spectrum as an inverse problem.We demonstrate the feasibility of the method in a well-defined model environment and investigate its opportunities and limitations for different types of input tracers.
The article is organized as follows.Section 2 recalls the fundamentals of the theory behind the age spectrum, makes explicit its relation to tracer concentrations and reviews previous approaches used to infer the age spectrum from tracers.Then, Section 3 presents the proposed inversion methodology and evaluates it based on idealized age spectra and a set of artificial decaying tracers.In Sect.4, the method is used on realistic age spectra from a chemistry transport model, which motivates a discussion of its limitations and applicability to observable tracers.Finally, Section 5 provides the conclusions.
2 Theoretical background: Relationship between age spectrum and tracers

Lagrangian path distribution
In the Lagrangian view of atmospheric transport (large-scale advection and mixing), each air parcel can be conceptually decomposed into an infinitude of infinitesimal and irreducible "fluid elements" that maintain their integrity against mixing for all timescales (Waugh and Hall, 2002).To each "fluid element" corresponds one Lagrangian path connecting a source (located at a given position on a surface) and the air parcel.Note that for any given source and emission time there might be a number of Lagrangian paths and hence of fluid elements.Each fluid element then explains a fraction m k of the mass of the air parcel, so that the partition of fluid elements fulfills: (1) Such a decomposition enables to understand the properties of the air parcel by disentangling the relative contribution of air masses of different origins.For instance, let us consider the age of the air parcel (average transit time since leaving the surface Ω) τ .This age of air can be broken down into the transit times of each of the fluid elements τ k , with the relation: Similarly, for a tracer of mixing ratio χ, one formally may write: It should be noted here that the decomposition used in Eq. 3 is not meaningful for all tracers.Actually, Eq. 3 makes sense only if the evolution of χ due to chemistry (or any process other than transport) can also be decomposed as: where is the rate of change of χ k within each of the k fluid elements, if they were separated (unmixed).In other words, dχ dt does not depend on whether the fluid elements are mixed or remain isolated from one another.For instance, reactive chemical species involved in 2 molecules reactions do not meet that requirement because their rate of change is, in general, affected by mixing (if the different fluid elements have different tracer concentrations).A simple example of tracers fulfilling the condition expressed in Eq. 4 is conserved tracers, for which dχ dt = dχ k dt k = 0. Another example is that of tracers with a linear decay or growth rate in their concentration (Schoeberl et al., 2000), such as radioactive tracers or tracers subject to photochemical loss.Their mixing ratio verifies: with r and t representing an eventual dependency of the growth/decay coefficient λ on position and time (for photochemical loss).For a pool of n tracers the description in Eq. 5 can easily be generalized to with χ the vector of trace species' mixing ratios and M(r, t) the matrix of growth/decay coefficients.Besides the "parent" radioactive tracers of Eq. 5, Equation 6 also encompasses the products of their decay ("daughter" tracers).
Being the probability distribution of transit times τ k for all fluid elements constitutive of the air parcel, the age spectrum can be viewed as a specific regrouping of Lagrangian paths according to transit time.It is also a boundary propagator of the continuity equation of conserved tracers from the surface Ω into the atmosphere (e.g., Holzer and Hall, 2000): in other words, the age spectrum relates the concentration of an inert tracer within the atmosphere to its uniform boundary condition on Ω.
This result can be extended to include tracers undergoing radioactive or chemical loss, as shown by a number of studies (e.g.Schoeberl et al., 2000;Waugh et al., 2003;Schoeberl et al., 2005).In the next subsection, we recall those analytical relations between age spectra and tracer content.This formal description will allow the reader to clearly apprehend the suitability of given trace gas species to probe the age spectrum.
2.2 Relation between age spectrum and tracer content 2.2.1 From age spectrum to tracer content: the forward model Assuming a uniform boundary condition in the surface region Ω, the mixing ratio ξ of a tracer with decay rate λ may be expressed as (e.g., Waugh et al., 2003): where τ is the transit time from Ω to (r, t) and ξ Ω is the tracer concentration at the surface.Here, G(τ ; r, t) represents the distribution of transit times, i.e. the age spectrum, at time t and position r.In the following we will drop the explicit reference to r in order to simplify the notations.Equation 7 can be generalized to a vector equation for n different tracers (similar to our arguing regarding Eq. 6): where bold fonts refer to vectors, M is the matrix of growth/decay coefficients and e Mτ the matrix exponential of Mτ .Equation 8 may encompass parent radioactive tracers as well as the whole associated decay chain (primary, secondary,... decay products), as explained in more detail in Sect.4.2.
Equations 7 and 8 show that the mixing ratios of conserved and linearly decaying (or growing) tracers may be expressed as the convolution of a generic function (involving time dependency of the source and chemistry) and the age spectrum G.
The tracer content can hence be seen as a weighted average of the age spectrum, and the functions e Mτ ξ Ω (t − τ ) as weighting functions (note that this perspective differs from the general view of the age spectrum as a weighting function of the tracer boundary condition history).However, although information on the age spectrum is contained in the tracer concentrations, it is far from being directly accessible because of this convolution with the tracer-dependent weighting functions.This limitation is evident in the case of linearly decaying tracers with constant boundary condition at the surface and constant lifetimes For those, Equation 8 simplifies as: where, as noted by Schoeberl et al. (2000), G is the Laplace transform of the age spectrum.(Note that with the normalization of ξ k (t) by its time-dependent surface value ξ Ω k (t), Equation 9 also applies to inert tracers exponentially increasing at the surface with growth rates 1/τ k .)The corresponding exponential weighting functions are represented in Fig. 1 for a pool of such tracers.They all peak for short transit times, so that the information provided by the different tracers is partly redundant and needs to be deconvolved.In general, this deconvolution may be achieved using different approaches, which will depend on the type of tracer considered and its associated weighting functions.

Diagnosing the age spectrum from the tracers: review of previous approaches
There have been a few attempts to characterize the age spectrum from the knowledge of tracer concentrations.Andrews et al. (1999) used time series of CO 2 and N 2 O to diagnose the transit time distribution, assumed to be a superposition of two inverse Gaussians.Johnson et al. (1999) used water vapor time series from which they deconvolved the age spectrum by the mean of Fourier transform.However, both studies heavily relied on the assumed stationarity of the atmospheric flow.In general, stratospheric transport and the associated stratospheric age spectrum are non-stationary, as evident from model simulations (e.g.Li et al., 2012;Diallo et al., 2012;Ray et al., 2014;Ploeger and Birner, 2016).In particular, the age spectrum exhibits seasonal and interannual variability.A few techniques have been proposed to estimate the age spectrum from tracer mixing ratios which do not rely on the stationarity assumption.They are briefly reviewed in the following.
A first approach, which might be referred to as the moment-estimate approach, is exposed for instance in Waugh et al. (2003).It is based on the relation between the moments of the age spectrum M n : and the concentration χ of a passive tracer (λ = 0 in Eq. 7) with a boundary condition χ Ω (t) evolving as a polynomial function of time t, of order N , so that one may write: where the α n are the coefficients of the polynomial.The relation is: Linearly increasing tracers constitute a particular case of Eq. 11 with N = 1 and α 1 < 0, such that χ Ω (t − τ ) = χ Ω (t) + α 1 τ .
For those, Equation 12 implies that the delay time ∆τ = (χ(t)−χ Ω (t))/α 1 is also the first moment of the age spectrum, called the mean age (e.g.Waugh and Hall, 2002) This last relation has been extensively used to derive mean age of air from linearly increasing conserved tracers, such as SF 6 and CO 2 (e.g.Engel et al., 2009).More generally, the moment-estimate approach builds on Eq. 12 to constrain specific moments of the age spectrum from tracers with different time dependency (linear, quadratic, ...).Knowledge of given moments (e.g. the first 2 moments) then enables to characterize the full age spectrum, assuming that the distribution has a given shape (such as an inverse Gaussian, solution of the age spectrum for 1D advection diffusion problems, as was done by Hall et al., 2002).This reasoning is however limited by the fact that real age spectra may exhibit a variety of shapes, and are not necessarily inverse Gaussians.
A second approach is the Boundary Impulse Response (BIR) method (Li et al., 2012).This method is based on a set of conserved pulse tracers, i.e. tracers which satisfy the boundary condition: In that case, the relation between the pulse tracer mixing ratio and the age spectrum reads Thus, a set of N such tracers initialized following Eq.14 with different, regular time intervals (i.e.t k = kδt) provides a resolved (though discretized) description of the age spectrum for transit times up to N δt.The BIR method has recently been employed in atmospheric chemistry-transport models (Li et al., 2012;Ploeger and Birner, 2016) in order to gain knowledge on the model age spectrum.Though a useful diagnostic in models, the BIR method requires this specific pool of artificial pulse tracers and cannot in general be applied to standard tracers that might be available from observations.
A third approach consists in retrieving (e.g. through least-square regression) the parameters of a given function representing the age spectrum so that it best fits the tracer concentrations.We will call that approach the parametric approach (Hall et al., 2002).Like the moment-estimate approach, it is based on the expectation that G has a given shape (e.g. an inverse Gaussian Hall and Plumb, 1994).The technique can easily be applied to observed tracers and was employed by Schoeberl et al. (2005).Although it provides reasonable results, the parametric approach suffers from the same caveat mentioned above that the shape needs to be assumed a priori.Very recent results show that it can be substantially improved by including information about the seasonality in transport (Hauck, M., Fritsch, F., Garny, H. and Engel, A.: Deriving stratospheric age of air spectra using chemically active trace gases, in prep.for Atmospheric Chemistry and Physics).
3 Inversion of the age spectrum from (non-pulse) tracers As emphasized by the review of the literature in the previous section, retrieving the age spectrum without assuming either stationarity of the flow or an a priori shape has never been attempted to our knowledge, although it has been suggested by some authors, including Schoeberl et al. (2000).Below, we describe a methodology to perform such retrievals and investigate its relevance for estimating the age spectrum.

Formulation of the discretized problem
Following Schoeberl et al. (2000), we discretize the convolution integral in Eq. 8 in transit time intervals [t j , t j+1 ] with the k-subscript indicating the k-th component of the tracer species vector and the "weighting function matrix" elements L k j and age spectrum vector G j given by To obtain the second equality in Eq. 16, we have assumed that G is piecewise constant over [t j , t j+1 ].We have also truncated the transit time axis at some t n , for practical computation reasons.In order to clarify the notation, we drop the explicit reference to t in G in the remainder of the paper, but it is implicit that the age spectrum depends on both time and location.Considering the full vector of mixing ratios, ξ, Eq. 16 can be written in matrix form For the special case of a suite of linearly decaying (radioactive) tracers with unit mixing ratio at the surface, as described by Eq. 9, the elements of the weighting function matrix L are simply A piecewise constant representation of the weighting function e M τ ξ Ω (t − τ ) could also have been used if no analytical expression had been available.
In order to gain information on G from the radioactive tracers, Schoeberl et al. (2000) suggested to use Eq.18 and to construct a square matrix L from which one could estimate G as G est = L −1 ξ obs .This method is not applicable in practice, however, because the problem is ill-posed and sensitive to small perturbation of ξ obs and because the matrix L is nearly singular (as demonstrated in appendix A1).

Inversion approach
Rather than directly inverting L, it is more appropriate to consider the determination of G from the observed trace gas mixing ratios ξ obs as an inverse problem, in which Eq. 18 is the forward model.In this formulation, the tracer content provides information on the convolution of the age spectrum with given functions.In that respect, it is similar to atmospheric soundings, for which the radiances measured at different wavelengths provide information on temperature and tracer profiles.Appropriate approaches to deal with such inverse problems are described in textbooks such as Rodgers (2000).In the following, we summarize the relevant pieces of information for the specific case considered here.
A solution to the discretized problem may be obtained through the minimization of a cost function J(G), here expressed as: The first term S −1 is the inverse covariance matrix of the "observed" (or modeled) tracers.It quantifies the departure from observations may correspond to instrumental noise or model error as well as uncertainties in the estimation of L (as, e.g., uncertainties in the decay coefficients, in the boundary condition ξ Ω k or even numerical errors).In our context, the second term involving the a priori G a and its inverse error covariance matrix S −1 a is introduced for regularization purposes (to avoid unphysical oscillations and large negative values of the retrieved G), in order to penalize solutions far from the a priori value.
Since the problem is already linear, the optimal G = G est (which minimizes J) can be readily estimated as: Contrary to most practical inverse problems, ours is of sufficiently small dimension (100 tracers and a few hundred points along the transit time axis at the most) so that a direct inversion of the matrix may be attempted without running into computational and memory limitations.However, similarly to most inverse problems, it is not obvious how to obtain values for the matrices S (which represents different sources of errors) and S a (which may only be estimated from models) nor to get a value for G a .
We will follow an empirical approach here, often referred to as Tikhonov regularization.We set G a = 0 and take S as σ 2 I and S a as σ 2 a α 2 I where I is the identity matrix, σ 2 a rough estimate of the variance of the "observation" (or model) error , σ 2 a a rough estimate of the variance of G a and α a positive scalar.Then the cost function can be rewritten: and the optimal estimate is: In practice, different values of α can be tested until a reasonable retrieval is obtained.For a certain range of α values, the retrievals are similar.That range encompasses the ratio of noise variance of the observation to the one of the a priori.
It should be noted that there is no guarantee that the estimated age spectra G est are positive for all transit times.As they are a result of optimal estimation, negative values should not be discarded, but taken into account in order to get the most accurate average and reduce the bias.

Feasibility and performance of the inversion
In order to test the feasibility of retrieving age spectra from a set of tracers, the sensitivity of the retrieval to noise in particular, preliminary checks with known, idealized spectra should be performed.We propose in this subsection a standard procedure to ensure the feasibility of the retrieval for a given tracer set and apply it to the particular case of the set of radioactive tracers presented in Fig. 1.

Idealized age spectrum and tracer set
The first step is to construct an age spectrum and the associated tracer composition as a test bed for the retrieval method.It is straightforward to estimate the decaying tracers from the perfect knowledge of the age spectrum, either analytically or through numerical integration of Eq. 18 with a fine resolution along the transit time axis.For the idealized age spectrum, we use the canonical expression for 1-D advective-diffusive systems given by (e.g.Waugh and Hall, 2002): where Γ is the mean age and ∆ the age spectrum width.This functional form for the age spectrum is known as inverse Gaussian function and has been extensively compared with model spectra (e.g.Schoeberl et al., 2005).The pseudo-observed (or modeled) mixing ratios of the tracers are derived as: Here, G hr j = G tj +tj+1 2 with t j = jδt and δt = 1 day.The error represents the uncertainty associated with the observation or modeling of the tracer.Since the accuracy of trace gas mixing ratios from measurements or models is generally proportional to their content, we take base proportional to the actual tracer mixing ratio, i.e. : where base is here a vector of random numbers from independent uniform distributions over [−0.05; 0.05].

Setting-up the retrieval
Typically, two parameters need to be chosen to set up a retrieval: the resolution along the transit time axis and the strength of the regularization, i.e. the value of α.It is actually advantageous to start with a high resolution along the transit time axis (e.g. 1 month) to nail down the value of α, before determining the effective resolution of the retrieval and adjusting the inversion to that resolution.
If the uncertainties associated with the observations or the a priori are not precisely known, there is some freedom in the choice of the optimal α.One procedure is to empirically test different values of α, and choose the best fit through visual inspection of the retrieved spectrum (i.e. until complete removal of the noise oscillations).However, this leaves room to subjectivity; a more objective approach is the L-curve optimality criterion (e.g.Hansen, 1992;Ungermann et al., 2011).This approach consists in plotting the residual obtained assuming different values of α in Eq. 23.For many inverse problems and a wide range of α, this yields a L-shaped curve, of which the corner (point of largest curvature) stands as a compromise between fidelity to the tracers and proximity to the a priori.
with G est calculated using Eq.23 with the corresponding value of α.The displayed curve is the average misfit-constraint for 100 retrievals from 100 sets of pseudo-observations with different realizations of the noise (different in Eq. 25).
In order to construct the L-shaped curve and to determine an appropriate value for α, we generate a set of 100 pseudoobservations ξ obs by varying in Eq. 25, with the "true spectrum" G hr taken as an inverse Gaussian with Γ = 2 years and ∆ = 1 year.For each of the 100 realizations of ξ obs , a retrieval G est is then performed using a given α in Eq. 23.This procedure is carried out for different values of α, resulting in 100 L-shaped curves (for each of the 100 realizations ξ obs ).
The average (for representativeness) of the resulting 100 L-shaped curves (i.e.turns out to be a good choice for our problem. Figure 3 shows the retrieved spectra obtained using α = 10 −1 , for two typical cases, a "young-age spectrum" (Γ = 2 years and ∆ = 1 year) and an "old-age spectrum" (Γ = 5 years and ∆ = 2 years).The thin black lines are individual retrieval results for 100 retrievals from the 100 sets of pseudo-observations including noise, while the thick black lines are the averages (shaded area: +/− 1 standard deviation).For both idealized spectra, the averages agree reasonably well with the input (red lines).In particular, the location of the mode is found in both cases and the general shape and magnitude of the spectrum are reproduced.
However, unrealistic negative values arise for small and large transit times (where the actual spectrum is close to 0), and the exact magnitude of the mode is not captured, with a ∼ 25% underestimation.Furthermore, there is a significant dispersion of individual retrievals around the average.This dispersion can be reduced by increasing the strength of the regularization α, but at the price of a deteriorated agreement of the multi retrieval average with the true spectrum.Conversely, a better agreement of the average spectrum with the input can be achieved by decreasing α, at the price of an increased dispersion in individual retrievals.As described above, the choice of α is a compromise between the reliability of individual retrievals and the accuracy of multi-retrieval averages.
We would like to emphasize that a different value of α may suit better when the relative strength of the noise is modified.
However, as the problem is ill-posed, regularization is required even in the absence of noise (see appendix A1).

Resolution
To perform the retrieval presented above, only 9 tracers were used whereas there were the 119 components of the spectrum to invert (monthly bins on a 10 year long transit time axis).It then comes without surprise that the retrieval is strongly underconstrained and requires regularization, especially since the weighting functions all peak at τ = 0.The effective resolution in transit time of the inverted spectrum can be investigated from the averaging kernel matrix A defined by Eq. 23 as The matrix A quantifies the contribution of the value of G at different transit times to the retrieved age spectrum G est at a specific transit time, and thus the resolution and ability to distinguish specific features.Averaging kernels peaking at one single transit time would provide the best resolution.For our setup, the averaging kernels are displayed in Fig. 4. As expected from the shape of the weighting function (Fig. 1), the resolution is better for short transit times, although even for those the effective resolution is worse than the 1-month-transittime bin size chosen for the retrieval, as can be seen from the overlap of the averaging kernels.The averaging kernels also exhibit negative lobes, which are responsible for the negative values seen in the retrieval at transit times characterized with low values of G.The amplitude of the negative values may be decreased by strengthening the regularization, but this reduces the sharpness of the peak of the averaging kernels and hence degrades the resolution.
Given the redundancy visible in the averaging kernels, it is possible to use a sparser resolution grid in transit time, which would better reflect the information available from the tracers.Although there is some freedom in the choice of the grid, we  McKenna et al., 2002;Konopka et al., 2004).The general setup of the model is described by Pommrich et al. (2014).The model simulation was started on 01/01/1979 and includes a pulse-tracer set to estimate the age spectrum using the BIR method similar to the one used by Ploeger and Birner (2016).From the pulse tracer mixing ratios the "true" model age spectra have been calculated independently using the BIR method, to validate the new age spectrum retrieval.In addition to the pulse tracers, 28 artificial radioactive tracers with boundary conditions at the surface and linear decay in the free atmosphere have been introduced.They consist in: one tracer with a decay time of 15 days, 17 with decay times ranging from 30 to 510 days with a 30-day step, and 10 with decay times from 570 to 1380 days with 90-day step.In Fig. 5, the age spectra retrieved from these linearly decaying tracers using the new method introduced in Sect. 3 are compared to age spectra estimated with the BIR method (Ploeger and Birner, 2016) for different altitude-latitude ranges on the 31/12/1983.
Figure 5 illustrates the unequal performance of the inversion in different cases.For short transit times, seen in the tropical upper troposphere (Fig. 5 a), the shape of the age spectrum is very well captured, despite the sharpness of the modal peak.At higher altitudes in the tropical pipe (Fig. 5 b), the transit time distribution exhibits two peaks, with the first mode corresponding to the (most recent) "direct ascension" from the surface while the second is a remainder from the increased entry of air in the stratosphere during the previous winter compared to the subsequent spring (see Ploeger and Birner, 2016, for further discussion of age spectrum seasonality).This bimodal behavior is smoothed out in the retrieval so that the two modes cannot be distinguished from one another in the retrieved spectrum, but the tail and general shape of the spectrum are well represented.
Only the mode from the previous winter has reached higher up (Fig. 5 c), which results in a translated spectrum with a larger tail compared to the ones displayed in panels a) and b).The full magnitude of the main peak is not reproduced in the retrieval, although its location is correct.The multipeak structure resulting from the seasonal cycle in the northern hemisphere stratosphere is completely smoothed out (Fig. 5 d).
The different examples above show that the (radioactively) decaying-tracer setup effectively enables to retrieve the general shape of the age spectrum.However, high-resolution features, such as the magnitude of individual peaks or the seasonal cycle in the age spectrum, are either underestimated or not retrieved at all, in particular fine-scale structures at large transit times.The comparison of the quality of the retrievals for different input spectra in panels a) and c) emphasizes the better resolution for short transit times.This is an immediate consequence of the shape of the averaging kernels, which are wider for large transit times (Fig. 4), due to the shape of the weighting functions for the radioactive tracers (Fig. 1).The resolution along the transit time axis is 1 month and the chosen regularization strength is The better quality of the retrievals for short transit times makes them most useful in the "ventilated" regions, i.e. the tropical pipe and the mid-latitude surf-zone.This is illustrated in Fig. 6, which contrasts the actual modal age of air determined with the BIR method with that derived using the retrieval procedure within the tropical pipe.As shown by Ploeger and Birner (2016), in the tropical pipe (as well as in the wintertime stratospheric "surf zone") the modal age is a useful indicator of the residual circulation transit time.Figure 6 shows that the retrieved tropical modal age agrees reasonably well with the BIR modal age (consistent with Fig. 5a, b).This is also the case for the "young" age spectra of the Midlatitude lower stratosphere (Fig. 5 c).However, for age spectra with long tails towards large transit times and a number of distinct peaks corresponding to the seasonal cycle, such as encountered in the midlatitude polar mid-stratosphere (Fig. 5 d), large errors occur.These errors partly originate from the coarser description of the spectrum at large transit times and partly from the inability of the inversion to capture the annual cycle.Introducing other tracers in the retrieval might allow to improve on this aspect, as investigated in the

Use of additional tracers to retrieve realistic age spectra
Besides parent radioactive tracers, Eq. 8 also encompasses daughter radioactive tracers, which are for instance the products of the decay of a surface emitted tracer, following the decay chain: The rate of change of the mixing ratio ξ B of the daughter tracer is given by: Let us now consider a set of parent and daughter tracers, with ξ A and ξ B the vectors of their mixing ratios.If the boundary condition at the surface for the parent tracers is ξ A = ξ Ω A and for the daughter tracers ξ Ω B = 0, and if the decay times are equal for each couple (τ Ak = τ B k = τ k ), then for each k the tracer mixing ratio is given by the following relation: We added a set of such daughter tracers to the parent radioactive tracers used in Sect.3.2.The method employed to initialize the tracers and set up the retrievals is the same as the one presented in Sect.3.2, except that the basic spectrum is now given by: where ω = 2π year −1 is the angular frequency, A the amplitude and φ the phase of the annual cycle, and dτ is a normalization constant.This functional form, introduced by Hauck et al. (in prep. for ACP), is an adjustment of Eq. 24 allowing the inclusion of the seasonal cycle.
Figure 8 shows the results of the retrieval experiment.The input spectra (red curves) bear resemblance with the realistic spectra in Fig. 5 (c, d).In particular, they exhibit a clear annual cycle, evident from the annually repeating peaks.The default retrieval using only parent tracers (black curves) does not fit this pattern, and has essentially the same shape as for an input without seasonal variability (as in Fig. 3).With both parent and daughter tracers (green curves), the fit to the input spectrum is improved.In particular, the uncertainty is clearly reduced.However, the seasonal variability is still absent from the retrieval.
To retrieve the seasonal variability in the age spectrum, we include another type of tracers.These are pairs of conserved tracers subject to periodic boundary conditions, such as sinusoidal tracers varying as: where ω m is the angular frequency of the oscillations.Such a pair of tracers in phase quadrature will provide additional information on periodic variations in the spectrum, like the seasonal cycle.The transfer matrix coefficients for those tracers are (calculated from Eq. 17) We added a set of sinusoidal tracers with periods of one and two years in addition to the set of parent radioactive tracers used in Sect.3.2 and the daughter tracers discussed above, to further improve the retrieval.The retrieval results are shown in Fig. 9.The addition of the periodic tracers (red curve) enables to retrieve the seasonality in the spectrum without deteriorating the representation of the general shape of the spectrum.Hence, it appears that with an adequate pool of time-varying tracers, high frequency features in the spectrum, such as the seasonal cycle, can be retrieved.

Application to observable tracers 10
It is beyond the scope of our study to retrieve atmospheric age spectra from actual tracer measurements.However, with appropriate considerations of errors in the definition of the error covariance matrix S (Eq.23), the methodology introduced in this paper can be straightforwardly applied to in situ or remotely observed tracer measurements.Tests with idealized distributions, as shown in Sect.3.2, enable to find out which properties of the transit time distribution can be inferred from a particular set of tracers.The experiments presented above already provide some general insight on this problem.Short-lived species with exponential decay or conserved tracers increasing exponentially at the surface can give detailed information on the transit time distribution for rapid transport, as suggested by Fig. 5.They might be sufficient to retrieve the age spectrum in the free troposphere resulting from convective transport from the boundary layer.For the stratosphere, with longer transport time scales involved and delayed arrivals of air masses, the parent radioactive tracers still carry some information on the transit time distribution, but their usefulness is more limited.Especially, they alone cannot be used to retrieve the annual cycle in age of air.However, they might be combined with long-lived tracers that exhibit an annual cycle (such as CO 2 ) to better constrain the age spectrum.

Conclusions 10
The concentrations in chemical tracers with different dependencies on transit time carry information on the age of air spectrum, the transit time distribution from the surface to a given location in the atmosphere.In this paper, we propose a method to retrieve the age of air spectrum from different trace gas species' mixing ratios.Formulating the question as an inverse problem, its dimension and complexity are by far more manageable than that of the inversions routinely performed for satellite retrievals of temperature and tracer profiles.In particular, the forward model (a mere convolution) is linear and the uncertainties are fairly well known compared to that of radiative transfer.Therefore, simple Tikhonov regularization appears sufficient to constrain the problem and retrieve the atmospheric transit time distribution.
Using prescribed age of air spectra and a set of artificial decaying radioactive tracers, we demonstrated the feasibility of the approach: even in the presence of forward model uncertainties and noise, the retrieved distributions are in reasonable agreement with the input age spectra.Furthermore, we applied the method to atmospheric transport simulations with the reanalysis-driven CLaMS model; the age spectra retrieved from a set of parent decaying tracers compared relatively well with spectra derived using the Boundary Impulse Response method, especially regarding the general shape of the distribution.However, fine-scale features, such as the seasonal cycle in transit-time probability, could not be captured with only-decaying tracers due to the large width of the averaging kernels.We show that the caveat may be circumvented by including trace gas species with seasonally varying concentrations at the surface and daughter decaying species in the retrieval.
The methodology introduced in this work can be applied in a number of situations.First, it might prove useful for the estimation of age spectra in models.Indeed, the most commonly used method, the Boundary Impulse Response method (Li et al., 2012;Ploeger and Birner, 2016), requires an increasing number of tracers with increasing maximum resolved transit time, which is cumbersome and computationally expensive, especially in Eulerian models (Li et al., 2012).It has in particular the disadvantage of a constant resolution as a function of transit time, which leads to unnecessarily high resolution to describe the tail at long transit times.With a refined set of artificial tracers (combining pulse and non pulse tracers), the inversion approach may enable an accurate and resolved description of the age spectrum at a reasonable computational cost.
However, the age spectrum retrieval approach might be most useful when trying to estimate transit time distributions from observations.An important number of tracers with different lifetimes and surface tendencies can nowadays be measured by Aircore hanging below balloons (Membrive et al., 2017;Engel et al., 2017) and whole air samplers onboard aircraft (as was done in some recent campaigns, e.g.Pan et al., 2017;Jensen et al., 2017).Although more limited in resolution, some remote sensing instruments, such as GLORIA (Riese et al., 2014) or ACE-FTS (Bernath et al., 2005), can also retrieve an important number of relevant atmospheric species on which this approach could be applied.We hope that our methodology will pave the way for a more precise and global characterization of the age spectrum from observations in the future.   .Age spectrum from a direct inversion (Eq.A1) using a set of radioactive tracers (black) vs input age spectrum (red).Note that due to the huge amplitude of the characteristic oscillations associated with the ill-posed, underconstrained problem, different y-axes are used for the inversed and input age spectra.
Author contributions.TEXT

Figure 2 .
Figure 2. Average L-curve for the 1 month-resolution setup.Each point of this curve corresponds to a pair < 1 σ 2 (L G est − ξ obs ) T (L G est − ξ obs ) > vs < 1

Figure 3 .
Figure 3. Input age spectra (red) and average retrieved age spectra (black), for an input idealized age spectrum with Γ = 2 years, ∆ = 1 year (left) and Γ = 5 years, ∆ = 2 years (right).The average retrieved age spectra are averages of 100 retrievals from 100 set of pseudo-tracer observations χ obs (i.e. 100 different realizations of the noise in Eq. 25).The gray shading corresponds to +/-the standard deviation of the hundred retrievals and shows the dispersion-noise induced uncertainty.The thin grey curves are individual retrievals.

Figure 4 .
Figure 4. Averaging kernels to the age spectrum for different retrieved transit times with the high-resolution (1 month) retrieval.(Left) actual averaging kernels.(Right) Averaging kernels normalized by their respective maximum value.Note that the averaging kernel at a particular transit time is the respective row of the averaging kernel matrix.
Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-990Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 5 October 2018 c Author(s) 2018.CC BY 4.0 License.keep the linear grid spacing in the following because of its simplicity and the demonstrated feasibility of the retrievals in that setup.4 Applications and discussion 4.1 Application to model data A first application of the inversion method is to retrieve age spectra from tracers in model simulations.To demonstrate this, we use a transport simulation performed with the 3D version of the Chemical Lagrangian Model of the Stratosphere (CLaMS

Figure 5 .
Figure 5. Age spectra retrieved from artificial decaying radioactive tracers (black) versus spectra estimated using the BIR method (red).The spectra on the different panels correspond to the same CLaMS model simulation on 31 December 1983, in different altitude-latitude regions.

Figure 6 .
Figure 6.Profile of tropical (15 • S-15 • N) modal age in CLaMS simulation, determined using the BIR method (red) or retrieved using the procedure highlighted in Sect. 3 (black).The red shading corresponds to +/-the standard deviation of the mode of the BIR spectrum in the 15 • S-15 • N region, and is introduced to guide the eye regarding the range of variability.

Figure 7 .
Figure 7. Shape of the weighting function to the age spectrum t τ k e − t τ k for radioactive or chemical product tracers with lifetimes τ k equal to that of the parent species.The vertical lines show the location of the maxima of the weighting functions, reached at transit times τ = τ k .

Figure 8 .
Figure8.(Red) Input age spectra, defined using Eq.32 with A = 0.3, φ = π 2 (for left and right panels), Γ = 2 years and ∆ = 1 year (left panel) and Γ = 5 years, ∆ = 2 years (right panel).Retrieved age spectra using (black) parent decaying tracers only or (green) both parent and daughter tracers.The full lines are average retrieved age spectra over 100 retrievals from 100 set of pseudo-tracer observations.The grey and green shadings corresponds to +/-the standard deviation of the hundred retrievals and show the dispersion-noise induced uncertainty.

Figure 9 .
Figure9.(Red) Input age spectra, defined using Eq.32 with A = 0.3, φ = π 2 (for left and right panels), Γ = 2 years and ∆ = 1 year (left panel) and Γ = 5 years, ∆ = 2 years (right panel).Retrieved age spectra using (black) parent decaying tracers only or (red) daughter and parent tracers and two sets of periodic tracers with periods of 1 and 2 years.The full lines are average retrieved age spectra over 100 retrievals from 100 set of pseudo-tracer observations.The grey and green shadings corresponds to +/-the standard deviation of the hundred retrievals and show the dispersion-noise induced uncertainty.
Atmos.Chem.Phys.Discuss., https://doi.org/10.5194/acp-2018-990Manuscript under review for journal Atmos.Chem.Phys.Discussion started: 5 October 2018 c Author(s) 2018.CC BY 4.0 License.The spectrum estimated using that approach is shown in Fig.A1.It exhibits large oscillations associated with the ill-posed, underconstrained problem.These oscillations are also present for a regression (Eq.23 with α 2 = 0) without the regularization terms (not shown), demonstrating the necessity of the regularization.

Figure A1
Figure A1.Age spectrum from a direct inversion (Eq.A1) using a set of radioactive tracers (black) vs input age spectrum (red).Note that