Simulating age of air and the distribution of SF₆ in the stratosphere with the SILAM model

The paper presents a comparative study of age of air (AoA) derived from several approaches: a widely used passive-tracer accumulation method, the SF₆ accumulation, and a direct calculation of an ideal-age tracer. The simulations were performed with the Eulerian chemistry transport model SILAM driven with the ERA-Interim reanalysis for 1980–2018.

The Eulerian environment allowed for simultaneous application of several approaches within the same simulation and interpretation of the obtained differences. A series of sensitivity simulations revealed the role of the vertical profile of turbulent diffusion in the stratosphere, destruction of SF₆ in the mesosphere, and the effect of gravitational separation of gases with strongly different molar masses.

The simulations reproduced well the main features of the SF₆ distribution in the atmosphere observed by the MIPAS (Michelson Interferometer for Passive Atmospheric Sounding) satellite instrument. It was shown that the apparent very old air in the upper stratosphere derived from the SF₆ profile observations is a result of destruction and gravitational separation of this gas in the upper stratosphere and the mesosphere. These processes make the apparent SF₆ AoA in the stratosphere several years older than the ideal-age AoA, which, according to our calculations, does not exceed 6–6.5 years. The destruction of SF₆ and the varying rate of emission make SF₆ unsuitable for reliably deriving AoA or its trends. However, observations of SF₆ provide a very useful dataset for validation of the stratospheric circulation in a model with the properly implemented SF₆ loss.

1 Introduction

The age of air (AoA) is defined as the time spent by an air parcel in the stratosphere since its entry across the tropopause (Li and Waugh, 1999; Waugh and Hall, 2002). The distribution of the AoA is controlled by the global atmospheric circulations, primarily the Brewer–Dobson and polar circulations. In particular, the temporal variation of AoA has been used as an indicator of the long-term changes in the stratospheric circulation (Engel et al., 2009; Waugh, 2009). AoA has been extensively used for evaluation and comparison of general circulation and chemical transport models in the stratosphere (Waugh and Hall, 2002; Engel et al., 2009).

Simulations of the AoA as defined above have been performed with Lagrangian transport models. The trajectories are initiated with positions distributed in the stratosphere and integrated backwards in time until they cross the tropopause. The time elapsed since the initialization is attributed as age of air at the point of initialization. Moreover, the distribution of the ages of particles originating from some location can be used to get the age spectrum there. Until recently, Lagrangian simulations of AoA did not explicitly account for turbulent mixing in the stratosphere (Eluszkiewicz et al., 2000; Waugh and Hall, 2002; Diallo et al., 2012; Monge-Sanz et al., 2012). Accounting for mixing adds up to 2 years to the mean AoA in the tropical upper stratosphere (Garny et al., 2014). In Lagrangian models, the mixing can be simulated with random-walk of the particles (Garny et al., 2014) or by inter-parcel mixing (Plöger et al., 2015; Brinkop and Jöckel, 2019).

The Eulerian simulations of AoA can be formulated in several ways. The approaches with an accumulating tracer, whose mixing ratio increases linearly in the troposphere, were used in a comprehensive study by Krol et al. (2018) and several studies before (e.g. Eluszkiewicz et al., 2000; Monge-Sanz et al., 2012). Another approach is to simulate a steady distribution of a decaying tracer, such as ²²¹Rn, emitted at the surface at a constant rate (Krol et al., 2018). Besides that, a special tracer that is analogous to the Lagrangian clock has been used. The tracer appears in the literature under names such as “clock-type tracer” (Monge-Sanz et al., 2012) or “ideal age” (Waugh and Hall, 2002). The ideal age has a constant rate of increasing of mixing ratio everywhere, except for the surface where it is continuously forced to zero. Similar tracers have long been used to simulate the transport times of oceanic water (e.g. England, 1995; Thiele and Sarmiento, 1990).

Direct observations of the age of air, as it is defined above, are not possible; therefore, AoA is usually derived from the observed mixing ratios of various tracers with known tropospheric mixing ratios and lifetimes (Bhandari et al., 1966; Koch and Rind, 1998; Jacob et al., 1997; Patra et al., 2011) or from the long-living tracers with known variations in the tropospheric mixing ratios. The studies published to date used carbon dioxide (CO₂; Andrews et al., 2001; Engel et al., 2009), nitrous oxide (N₂O; Boering et al., 1996; Andrews et al., 2001), sulfur hexafluoride (SF₆; Waugh, 2009; Stiller et al., 2012), methane (CH₄; Andrews et al., 2001; Remsberg, 2015), and various fluorocarbons (Leedham Elvidge et al., 2018).

For accumulating tracers, the mean AoA at some point in the stratosphere is calculated as a lag between the times when a certain mixing ratio is observed near the surface and at that point. The lag time is equivalent to the mean AoA defined above only in the case of the strictly linear growth and the uniform distribution of the tracer in the troposphere (Hall and Plumb, 1994).

In reality, there is no tracer whose mixing ratio in the troposphere grows strictly linearly. The violation of the assumption of the linear growth leads to biases in the resulting AoA distribution and its trends. It has been pointed out that the increasing growth rates of CO₂ and SF₆ lead to a low bias of AoA and its trends and make these tracers ambiguous proxies of the AoA (Garcia et al., 2011). Various corrections have been applied in several studies (Hall and Plumb, 1994; Waugh and Hall, 2002; Engel et al., 2009; Stiller et al., 2012; Leedham Elvidge et al., 2018) to deduce the “true” AoA from observations of tracers with the increasing growth rates. The effect of the correction method on the AoA estimates has not been investigated and must be considered a source of uncertainty in the resulting estimates. Thus, Garcia et al. (2011) concluded that accounting for the biases in the trend estimates due to varying growth rates would likely require uniform and continuous knowledge of the evolution of the trace species, which is not available from any existing observational dataset. Recently Leedham Elvidge et al. (2018) showed a minor sensitivity of the AoA to the choice of the correction method but without detailed analysis of the assumptions behind these methods. For a similar problem with the ages of oceanic water, it has been shown (Waugh et al., 2003) that, in the case of a inhomogeneously growing tracer, the tracer age is strongly influenced by the shape of the transient time distribution (TTD, also known as the “age spectrum”) at the particular location and time.

Another major source of uncertainty in the observational AoA is the violation of conservation of the tracer due to sources and sinks, such as oxidation of carbon monoxide and methane for CO₂ or mesospheric destruction for SF₆. The mesospheric sink of SF₆ leads to “over-ageing”, especially pronounced in the area of the polar vortices. The magnitude of the over-ageing was estimated to be as at least 2 years (Waugh and Hall, 2002). Besides being visible in many evaluations, e.g. Stiller et al. (2012, Fig. 4) and Kovács et al. (2017, Fig. 8), the over-ageing of the polar winter stratospheric air was studied by Ray et al. (2017, Fig. 4) within the dedicated exercise.

The simulations of SF₆ and the AoA in the atmosphere with the WACCM model (Kovács et al., 2017) have also reproduced the effect of over-ageing. However, its magnitude was much smaller than that inferred from the SF₆ retrievals of the limb-viewing MIPAS (Michelson Interferometer for Passive Atmospheric Sounding) instrument operated on board of the Envisat satellite in 2002–2012 (Stiller et al., 2012) and from the in situ observations of the ER-2 aircraft (Hall et al., 1999). Kovács et al. (2017) offered two possible reasons for the discrepancy: either SF₆ loss is still underestimated in WACCM or MIPAS SF₆ observations are low biased above ∼20 km. Neither of the cases have been analysed in depth, which leaves the status of MIPAS, currently the richest observational dataset for the stratospheric SF₆, unclear.

The aim of the present study is to provide self-consistent simulations of the spatio-temporal distribution of the AoA and of the SF₆ mixing ratio in the troposphere and the stratosphere during the last 39 years. The main modelling tool is the Eulerian chemistry transport model SILAM (System for Integrated modeLling of Atmospheric coMposition). The stratospheric balloon observations and retrievals of the limb-viewing MIPAS instrument mentioned above are used for validation of the simulated distribution.

With these simulations we

• compare different methods of estimating the AoA and quantify the inconsistencies in the AoA and its trends arising from violations of the underlying assumptions behind each method,

• analyse the causes of the discrepancies in the upper stratosphere between different methods of deriving the AoA,

• provide a solid basis for further studies of stratospheric circulation with observations of various trace gases and for studies of climate effects of SF₆.

The paper is organized as follows. Section 2 gives an overview of the modelling tools and the modelling and observational data used for the study. Section 3 describes the developments made for SILAM in order to perform the simulations: vertical eddy-diffusivity parameterization in the stratosphere and the lower mesosphere and the SF₆ destruction parametrization, as well as the model configuration used for the study. The sensitivity tests and evaluation of the simulations against the MIPAS retrievals and stratospheric balloon measurements of SF₆ mixing ratios are given in Sect. 4. Sensitivity of the AoA and its trends to the simulation setup and the choice of particular SF₆ tracer as an AoA proxy is studied in Sect. 5. The uncertainties of the used modelling approach and implications of AoA derived from SF₆ tracer are discussed in Sect. 6. The results are summarized in Sect. 7.

2 Methods and input data

2.1 SILAM model

SILAM (System for Integrated modeLling of Atmospheric coMposition, http://silam.fmi.fi, last access: 13 May 2020) is an offline 3D chemical transport model. SILAM features a mass-conservative positive-definite advection scheme that makes the model suitable for long-term runs (Sofiev et al., 2015). The model can be run at a range of resolutions starting from a kilometre scale in a limited-area up to a global coverage. The vertical structure of the modelling domain consists of stacked layers starting from the surface. The layers can be defined either in z- or hybrid sigma-pressure coordinates. The model can be driven with a variety of NWP (numerical weather prediction) or climate models.

The global 3D simulations of atmospheric transport for a variety of tracers representing AoA and SF₆ (see Sect. 3.4 for details) were performed with SILAM for the years 1980–2018 with the global long–lat grid of $\mathrm{1.44}{}^{\circ }\times \mathrm{1.44}{}^{\circ }$ cells (250×123 grid cells plus polar closures) and 60 hybrid sigma-pressure layers starting from the surface. The uppermost layer was between pressures of 0.1 and 0.2 hPa, whereas other layer bounds corresponded to the half levels of the meteorological driver – the ERA-Interim reanalysis (Sect. 2.2). The model time step was 15 min and the output consisted of daily-mean 3D concentrations of the tracers and air density. Emission data were taken from the SF₆ emission inventory (Rigby et al., 2010), which was extrapolated until 2016 as described in Sect. 3.4. Physical–chemical transformations of the SF₆-related tracers required developments described in Sect. 3.3.

In order to accurately model the AoA and the needed tracers, the vertical diffusion part of the transport scheme of SILAM has been refined to account for gravitational separation. In addition, several tracers with corresponding transformation routines have been implemented into the model. The SILAM configuration, used for the present study, is described in Sect. 3.4.

2.2 ECMWF ERA-Interim reanalysis

The ERA-Interim reanalysis of the European Centre for Medium-Range Weather Forecasts (ECMWF) had been used as a meteorological driver for our simulations. The dataset has T255 spectral resolution and covers the whole atmosphere with 60 hybrid sigma-pressure levels having the uppermost layer from 0.2 to 0 hPa with nominal pressure of 0.2 hPa (Dee et al., 2011). The reanalysis uses a 12 h data assimilation cycle, and the forecasts are stored with a 3 h time step. We used the fields retrieved from the ECMWF's MARS archive on a long–lat grid, 500×250 points, with a step of 0.72^∘. The four forecast times (+3, +6, +9 and +12 h) were used from every assimilation cycle to obtain a continuous dataset with 3 h time step. To drive the dispersion model, the data on horizontal winds, temperature, and humidity for 1980–2018 were used.

Since the resolution of the driving meteorology was twice higher than that of SILAM, the meteorological input for both cell interface for winds and cell mid-points for other parameters (surface pressure, temperature, and humidity) was available without interpolation. The gridded ERA-Interim fields are, however, a result of reprojection of the original meteorological fields computed as spherical harmonics. Moreover, the difference in the topmost layer of the ERA-Interim and SILAM data required vertical reprojection at the top of the domain. Together with the limited precision of the gridded fields retrieved from the ECMWF archive, they caused some inconsistency between the surface-pressure tendencies and the vertically integrated air-mass fluxes calculated from the meteorological fields in SILAM. Albeit small, such inconsistencies cause spurious variations in wind-field divergence that might result in gradual accumulation of errors in the tracer mixing ratios. To maintain strict global and local air-mass budget throughout the run, the wind fields were adjusted by distributing the residuals of pressure tendency and vertically integrated horizontal air-mass fluxes as a correction to the horizontal winds, as suggested by Heimann and Keeling (1989). The correction was, at most, of the order of centimetres per second, which is comparable to the precision of the input wind fields. The vertical wind component was then rediagnosed from the divergence of the horizontal air-mass fluxes for the SILAM layers as described in Sofiev et al. (2015). Validity of this procedure was demonstrated by its authors Heimann and Keeling (1989) and its applicability to the current case was confirmed in the Sect. 5.2 by comparison with another model simulations driven by ERA-Interim (Diallo et al., 2012).

The ERA-Interim reanalysis has been used earlier for Lagrangian simulations of AoA (Diallo et al., 2012) and has been found to provide ages that agree with those inferred from in situ observations in the lower stratosphere.

2.3 MIPAS observations of SF₆

To evaluate the results of the SF₆ modelling, we used the data from the MIPAS instrument operated on board Envisat in 2002–2012. MIPAS is a limb-sounding Fourier transform spectrometer with a high spectral resolution measuring in the infrared part of spectrum. Due to its limb geometry, the instrument provided good vertical resolution of the derived trace-gas profiles and showed high sensitivity to low-abundance species around the tangent point. Along the orbit path, MIPAS measured a profile of atmospheric radiances about every 400 km with an altitude coverage, in its nominal mode, from 6 to 70 km. The vertical sampling was 1.5 km in the lower part of the stratosphere (up to 32 km) and 3 km above, with a vertical field of view covering 3 km at the tangent point. Over a day, about 1300 profiles along 14.4 orbits were measured, covering all latitudes up to the poles at sunlit and dark conditions. The vertical distributions of trace gases were derived from the radiance profiles by an inversion procedure, fitting simulated spectra to the measured ones while varying the atmospheric state parameters.

The retrieval of SF₆ is based on the spectral signature of this species in the vicinity of 10.55 µm wavelength and is described in Stiller et al. (2008), Stiller et al. (2012), and Haenel et al. (2015). In the current study, we use an updated version of the SF₆ data (compared to the one described in Haenel et al., 2015) called V5H/R_SF6_21/224/225. The new algorithm uses the new absorption cross-section data on the SF₆ and a new CFC-11 band in the vicinity of the SF₆ signature by Harrison (2018) instead of the older cross-section data by Varanasi et al. (1994). The updated version provides up to 0.6 pmol mol⁻¹ higher SF₆ mixing ratios in the upper part of the stratosphere (above 30 km) than the old versions and is closer to independent reference data. Note that whilst we regard this newer version of MIPAS SF₆ data as an improvement, it has not yet been reported in a publication, and on that basis it is subject to uncertainty.

The retrieved profiles are sampled on an altitude grid spaced at 1 km, whereas the actual resolution of the profiles is between 4 and 10 km for altitudes below 30 km. The retrievals are supplemented with averaging kernels and error covariance matrices describing the uncertainties due to random noise in the radiance measurements, hereinafter referred to as measurement noise error, target noise error, or retrieval noise error. This error component, which is normally of the order of 10 % of the retrieved value, is fully uncorrelated from profile to profile, and therefore it virtually cancels out when averaged over a large number of profiles. In contrast, there exist systematic error components that are fully correlated between the profiles. Their assessment is difficult and depends on the knowledge about sources of systematic errors. Stiller et al. (2008) has assessed them to be of the order of 10 % at 60 km and 4 % at 30 km. These error components have to be considered when comparisons of monthly or seasonal means with other data are performed.

3 SILAM developments

Destruction of atmospheric SF₆ occurs at altitudes above 60 km (Totterdill et al., 2015) that fall within the topmost layer of the ERA-Interim data. The exchange processes in the upper stratosphere and lower mesosphere have to be adequately parameterized together with the destruction process. In our simulations we have suppressed the transport of SF₆ with mean wind through the modelling domain top (0.1 hPa, 65 km) and parameterized the SF₆ loss due to the eddy and molecular diffusion towards the altitudes where the destruction occurs. In this section we introduce the set of parameterizations that were implemented in SILAM for this study.

Figure 1 Vertical profiles of diffusion coefficients are shown. The distribution of the ERA5 profiles of the “mean turbulent diffusion coefficient for heat” parameter, molecular diffusivity for SF₆ in the U.S. Standard Atmosphere, and the three prescribed K_z profiles are shown. The eddy diffusion profile due to breaking gravity waves (after Lindzen, 1981) is given for the reference.

3.1 Eddy diffusivity

A large variety of vertical profiles for eddy diffusivity in the stratosphere and the lower mesosphere can be found in literature. In many studies in the 1970s–1980s, the vertical profiles were derived from observed tracer concentrations neglecting the mean transport. Most studies suggested that the vertical eddy diffusion has a minimum of 0.2–0.5 m² s⁻¹ (Pisso and Legras, 2008) at 15–20 km, agreeing quite well to the ones derived from the radar measurements in the range of 15–20 km (Wilson, 2004). Above that altitude, K_z was suggested to gradually increase by about 1.5 orders of magnitude towards 50 km due to breaking gravity waves (Lindzen, 1981).

The theoretical estimates of the effective exchange coefficients, considering the layered and patchy structure of stratospheric turbulence, suggest 0.5–2.5 m² s⁻¹ for the upper troposphere and 0.015–0.02 m² s⁻¹ for the lower stratosphere (Osman et al., 2016), which is about an order of magnitude lower than the estimates above.

The values of the eddy exchange coefficient at heights of 10–20 km estimated from the high-resolution balloon temperature measurements (Gavrilov et al., 2005) are ∼0.01 m² s⁻¹ with no noticeable vertical variation. It is not clear, however, how representative the derived values are for UTLS (upper troposphere and lower stratosphere) in general. We could not find any reliable observations of vertical diffusion in a range of 30–50 km.

The parameterization for vertical eddy diffusivity above the boundary layer used in SILAM has been adapted from the IFS model of the European Centre for Medium-Range Weather Forecasts (ECMWF, 2015). However, in the upper troposphere the predicted eddy diffusivity is nearly zero. For numerical reasons, a lower limit of 0.01 m² s⁻¹ is set for K_z in SILAM. Our sensitivity tests have shown that long-term simulations are insensitive to this limit as long as it is low enough. The K_z in the stratosphere is routinely set to the limiting value with relatively rare peaks, mostly in UTLS. Such a scheme essentially turns off turbulent diffusion in the stratosphere. The same is true for the recent ERA5 reanalysis dataset (Copernicus Climate Change Service (C3S), 2017) that provides the values of K_z among other model-level fields: the eddy diffusion routinely falls below the molecular diffusivity above 40 km (Fig. 1).

As a reference for this study, we took a tabulated profile of Hunten (1975), as it was quoted by Massie and Hunten (1981). The original profile covers the range up to 50 km, and the extrapolation up to 80 km matches the theoretical estimates by Lindzen (1981) and by Allen et al. (1981). We approximate the profile as a function of pressure in the range of 100–0.01 hPa (15–60 km):

$$\begin{array}{}\text{(1)}& K_z\left(p\right)=\mathrm{8}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}{\left({\displaystyle \frac{\mathrm{1}\mathrm{hPa}}{p}}\right)}^{\mathrm{0.75}}.\end{array}$$

The approximated profile was stitched with the default SILAM profile with a gradual transition within an altitude range of 10–15 km to keep the tropospheric dispersion intact. This profile gives values of K_z 3–6 orders of magnitude higher than the ones provided by the ERA5 reanalysis (Fig. 1) and 1–2 orders of magnitude higher than the estimates of Legras et al. (2005).

In order to cover the range of K_z values between the ERA5 profiles and the reference one (Eq. 1), we used two intermediate profiles obtained by scaling the reference one with factors 0.03 and 0.001. The three prescribed eddy-diffusivity profiles are hereinafter referred to as “1-Kz”, “0.03-Kz”, and “0.001-Kz”, respectively. The dynamic eddy-diffusivity profile adopted from the ECMWF IFS is referred to as “ECMWF-Kz”. In all simulations, the parameterization of K_z in the troposphere is the same, and linear transition from the SILAM K_z to the prescribed one occurs in the altitude range of 10–15 km.

3.2 Molecular diffusivity and gravitational separation

In tropospheric and stratospheric chemistry transport models (CTMs), gaseous admixtures are transported as tracers (i.e. advection and turbulent mixing do not depend on the species properties), whereas the molecular diffusion is negligible. Models that cover the mesosphere, such as WACCM (Smith et al., 2011), account for molecular diffusion explicitly. Since some of the K_z parameterizations of the previous section often result in values below the molecular diffusivity, the parametrization of molecular diffusion has been implemented in SILAM.

The molecular diffusivity of SF₆ in the air at temperature T₀=300 K and pressure p₀=1000 hPa is $D_{\mathrm{0}}={\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ (Marrero and Mason, 1972, Table 22). The diffusivity at different temperature T and pressure p is given by

$$\begin{array}{}\text{(2)}& D=D_{\mathrm{0}}{\displaystyle \frac{p_{\mathrm{0}}}{p}}{\left({\displaystyle \frac{T}{T_{\mathrm{0}}}}\right)}^{\mathrm{3}/\mathrm{2}},\end{array}$$

(e.g. Cussler, 1997). The vertical profile of molecular diffusivity in the U.S. Standard Atmosphere (NOAA et al., 1976) is shown in (Fig. 1). Note that the value for the reference diffusivity of SF₆ used in this paper is about a half of the one used in simulations with WACCM by Kovács et al. (2017). The reason is that WACCM uses a universal parametrization (Smith et al., 2011, Eq. 7 there) for all compounds. That parametrization relies solely on molecular mass of a tracer and does not account for, for example, the molecule collision radius. The latter is about twice larger for SF₆ than for most of stratospheric tracers. Thus, for this study we use the value from Marrero and Mason (1972), which results from fitting laboratory data for diffusion of SF₆ in the air.

The vertical diffusion transport velocity of admixture with number concentration $\stackrel{\mathrm{̃}}{n}$ and molecular mass $\stackrel{\mathrm{̃}}{\mathit{\mu }}$ in neutrally stratified media is given by Mange (1957):

$$\begin{array}{}\text{(3)}& w=-D\left[{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\mu }}}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\mu }}}{\partial z}}+\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\mu }}}{\mathit{\mu }}}-\mathrm{1}\right){\displaystyle \frac{\mathit{\mu }g}{kT}}\right],\end{array}$$

where μ is molecular mass of air, g is acceleration due to gravity, k is the Boltzmann constant, and T is temperature. With the ideal gas law p=nkT, in which p is pressure and n is number concentration, and the static law $\mathrm{d}p/\mathrm{d}z=-g\mathit{\rho }$, where ρ=μn is air density, Eq. (3) can be reformulated in terms of admixture mixing ratio $\mathit{\xi }=\stackrel{\mathrm{̃}}{n}/n$ and pressure. Then the vertical gradient of the equilibrium mixing ratio will be

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial \mathit{\xi }}{\partial p}}=\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\mu }}}{\mathit{\mu }}}-\mathrm{1}\right){\displaystyle \frac{\mathit{\xi }}{p}}.\end{array}$$

It is non-zero for an admixture of a molecular mass different from the one of air. Integrating the gradient Eq. (4) over the vertical, one can find that the equilibrium mixing ratios ξ₁ and ξ₂ at two levels with corresponding pressures p₁ and p₂ are related as

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{{\mathit{\xi }}_{\mathrm{1}}}{{\mathit{\xi }}_{\mathrm{2}}}}={\left({\displaystyle \frac{p_{\mathrm{1}}}{p_{\mathrm{2}}}}\right)}^{\stackrel{\mathrm{̃}}{\mathit{\mu }}/\mathit{\mu }-\mathrm{1}}.\end{array}$$

For heavy admixtures, such as SF₆ ($\stackrel{\mathrm{̃}}{\mathit{\mu }}=\mathrm{0.146}$ kg mol⁻¹) the equilibrium gradient of a mixing ratio is substantial. For example, the difference of the equilibrium mixing ratio of SF₆ between 0.1 and 0.2 hPa is a factor of 16.

In most of the atmosphere, the effect of gravitational separation is insignificant due to the overwhelming effect of other mixing mechanisms, whereas in the upper stratosphere the molecular diffusivity may become significant. Therefore, in the upper stratosphere heavy gases can no longer be considered tracers and the molecular diffusion should be treated explicitly. The effect of gravitational separation of nitrogen and oxygen isotopes in the stratosphere has been observed (Ishidoya et al., 2008, 2013; Sugawara et al., 2018); however, for isotopes the ratio of masses is relatively small, so the observed differences were also small (up to 10⁻⁵). For SF₆, the molecular mass difference is much larger.

In order to enable the gravitational separation in SILAM, we have introduced the molecular diffusion mechanism, which can be enabled along with the turbulent diffusion scheme. The exchange coefficients due to molecular diffusion between the model layers are precalculated according to Eq. (3) and discretized for the given layer structure for each species according to its diffusivity and molar mass. The U.S. Standard Atmosphere (NOAA et al., 1976) was assumed for the vertical profiles of temperature and air density during precalculation of the exchange coefficients. The exchange has been applied throughout the domain at every model time step with a simple explicit scheme.

3.3 SF₆ destruction

Estimates of AoA from the SF₆ tracer rely on the assumption of it being a passive tracer. SF₆ is indeed essentially stable in the troposphere and the stratosphere. IPCC (2013, Sect. 8.2.3.5) mentions that photolysis in the stratosphere as the main mechanism of SF₆ loss but without any reference to original studies. The statement is probably taken from Ravishankara et al. (1993). Reddmann et al. (2001) pointed at associative electron attachment in the upper stratosphere and mesosphere as the main destruction mechanism for SF₆ below 80 km. The recent study of Totterdill et al. (2015) gives some 1–2 orders of magnitude slower rates of electron attachment but keeps it the dominant mechanism of the SF₆ destruction in the altitude range up to 100 km. The highest destruction rate of 10⁻⁵ s⁻¹ occurs at the altitude of 80 km (Fig. 2). An important feature of this profile is that the destruction rate becomes significant above the top of our modelling domain (0.1 hPa, 65 km). The ERA-Interim meteorological fields have the uppermost level at 0.1 hPa and do not resolve the vertical structure of the atmosphere above that level. In order to assess the loss of SF₆, we have to parameterize the combined effect of the SF₆ transport through the 0.1 hPa and its destruction. Then the resulting fluxes can be applied as the upper boundary condition for our simulations.

Figure 2 The vertical profiles of SF₆ destruction rate (after Totterdill et al., 2015) and its approximation in the range of 55–75 km, given by Eq. (6).

As an approximation to the vertical profile of the destruction rate in an altitude range of 50–80 km, we have fitted the corresponding part of the curve in Fig. 9a of Totterdill et al. (2015) with a power function of pressure (magenta line in Fig. 2):

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{8}}{\mathrm{s}}^{-\mathrm{1}}{\left({\displaystyle \frac{\mathrm{0.2}\mathrm{hPa}}{p}}\right)}^{\mathrm{3}},\end{array}$$

where τ is the lifetime of SF₆ at the altitude corresponding to pressure p.

Figure 3 Vertical profiles of steady-state upward flux of SF₆ normalized with mass mixing ratio, F(p)∕ξ(p), for eddy-diffusivity and lifetime profiles given by Eqs. (1) and (6). The topmost model layer of SILAM and effective lifetimes of SF₆ there due to the destruction in the mesosphere for different K_z profiles are given.

The topmost level of the ERA-Interim meteorological dataset is located at 0.1 hPa, which is below the layer where the destruction of SF₆ occurs. Therefore, we have to put a boundary condition on our simulations to account for the upward flux of SF₆ through the upper boundary of the simulation domain. For that, we assume that the SF₆ distribution above the computational domain top is in equilibrium with the destruction and the vertical flux.

Assuming the profiles for K_z(p) and the SF₆ lifetime τ(p) are given by Eqs. (1) and (6), one can obtain a steady-state distribution of the mass-mixing ratio, ξ, of SF₆ due to destruction in the mesosphere at any point where both Eqs. (1) and (6) are valid and vertical advection is negligible. The latter assumption implies that the diffusive vertical flux overwhelms the advective one. The validity and implications of neglecting the regular vertical transport are discussed below. The steady-state profile of ξ can be obtained from a solution of the steady-state diffusion equation with a sink:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial \mathit{\xi }}{\partial t}}=g{\displaystyle \frac{\partial }{\partial p}}\left(F\right)-{\displaystyle \frac{\mathit{\xi }}{\mathit{\tau }\left(p\right)}}=\mathrm{0},\end{array}$$

where ρ(p) is air density, g is acceleration due to gravity, and the upward flux of SF₆ is given by

$$\begin{array}{}\text{(8)}& F\left(p\right)=g{\mathit{\rho }}^{\mathrm{2}}{K}_z\left(p\right){\displaystyle \frac{\partial \mathit{\xi }}{\partial p}}.\end{array}$$

The above equation was solved numerically as a boundary value problem with unit mixing ratio at a height of 1 hPa and vanishing flux, F(p) at p=0, for the set of K_z profiles. The shooting method with bisection was used to get the steady-state profiles of ξ(p) and F(p), corresponding to ξ(1 hPa)=1. For all considered cases, the flux F(p) decreased by several orders of magnitude already at the level of a few pascals (Pa), i.e. below the maximum of the depletion profile of Totterdill et al. (2015), indicating that the particular shape of τ(p) above that level does not influence the fluxes at the domain top (0.1 hPa). The steady-state upward flux of SF₆ F(p) normalized with the corresponding mixing ratio at each pressure, F(p)∕ξ(p), for the three test profiles of K_z is shown in Fig. 3 with solid lines.

The gravitational separation can be accounted for by introducing a term responsible for molecular diffusion and its equilibrium state Eq. (4) into the vertical flux Eq. (8) :

$$\begin{array}{}\text{(9)}& F\left(p\right)=g{\mathit{\rho }}^{\mathrm{2}}{K}_z\left(p\right){\displaystyle \frac{\partial \mathit{\xi }}{\partial p}}+g{\mathit{\rho }}^{\mathrm{2}}D\left(p\right)\left({\displaystyle \frac{\partial \mathit{\xi }}{\partial p}}-{\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathit{\mu }}-\mathit{\mu }}{\mathit{\mu }}}{\displaystyle \frac{\mathit{\xi }}{p}}\right).\end{array}$$

The profiles of F(p)∕ξ(p) resulting from F(p) in Eq. (7) are given in Fig. 3 with dashed lines. The magnitude of F(p)∕ξ(p) gives an equivalent regular vertical air-mass flux that would result in the same vertical flux of SF₆ if it were passive and non-diffusive. The equivalent regular vertical velocity ω_eq (in units of the Lagrangian tendency of a parcel pressure due to vertical advection) can be expressed as

$$\begin{array}{}\text{(10)}& {\mathit{\omega }}_{\mathrm{eq}}=-gF\left(p\right)/\mathit{\xi }\left(p\right).\end{array}$$

Accounting for molecular diffusion may either enhance or reduce the upward flux of SF₆ in the model. Along with setting the equilibrium state with the bulk of a heavy admixture being in the lower layers, molecular diffusion provides additional means for transport to the upper layers where the destruction occurs. For very low eddy diffusivities, the molecular diffusion is a sole mechanism of the upward transport of SF₆ towards depletion layers. For higher eddy diffusivity, the effect of molecular diffusion and gravitational separation becomes negligible.

For the model consisting of stacked well-mixed finite layers, the loss of SF₆ from the topmost layer due to the steady upward flux would be proportional to the SF₆ mixing ratio in the layer. This loss of mass is equivalent to a linear decay of SF₆ in the layer at a rate

$$\begin{array}{}\text{(11)}& {\mathit{\tau }}^{-\mathrm{1}}=g{\displaystyle \frac{F\left(p\right)}{\mathit{\xi }\left(p\right)\mathrm{\Delta }p}},\end{array}$$

where Δp is pressure drop in the layer.

In the upper layer of our simulations (between 0.1 and 0.2 hPa, grey rectangle in Fig. 3), the SF₆ lifetime τ due to turbulent diffusion is about 3 d for K_z of Eq. (1). After scaling the K_z(p) profile with factors of 0.03 and 0.001, one gets the lifetimes of 15 and 60 d, respectively. Note that the molecular diffusion sets the upper limit to the SF₆ lifetime in the topmost model layer: it can not be longer than 60 d for the 0.1–0.2 hPa layer. Close to this regime, the system becomes insensitive to the actual profile and values of the turbulent diffusion coefficient. The loss of SF₆ through the domain top was implemented as a linear decay of SF₆ in the topmost model layer, at a rate corresponding to the K_z(p) profile used in each simulation.

3.4 Simulated tracers

SILAM performs the 3D transport by means of a dimension split: transport along each dimension is performed separately as 1D transport. To minimize the inconsistency between the tracer transport and air-mass fluxes caused by the dimension split at finite time step, the splitting sequence has been inverted at each time step. The residual inconsistency was resolved by using a separate unity tracer, which was initialized to the constant mass mixing ratio of 1 at the beginning of a simulation. Should advection be perfect, the concentration of the unity tracer would be equivalent to air density (mixing ratio would stay equal to 1). The mixing ratios of the simulated tracers were then evaluated as a ratio of the tracer mass in a cell to the mass of the unity tracer.

In order to assess the effects of gravitational separation and destruction on the atmospheric distribution of SF₆, we used four tracers: SF₆ as a passive tracer sf6pass, SF₆ with gravitational separation but no destruction sf6nochem (no chemistry), SF₆ with destruction but no gravitational separation sf6nograv, and SF₆ with both gravitational separation and destruction in the upper model level sf6.

All SF₆ tracers had the same emission according to the SF₆ emission inventory (Rigby et al., 2010). The inventory covers 1970–2008 and was extrapolated with a linearly growing trend of 0.294 Gg yr⁻² until July 2016. The last 2.5 years were run without the SF₆ emissions to evaluate its destruction rate. Note that the emission extrapolation gives 9.4 Gg yr⁻¹ for 2016, which is somewhat higher than the later estimate of 8.8 Gg yr⁻¹ (Engel et al., 2018).

Besides the four SF₆ tracers, we used a passive tracer emitted uniformly at the surface at constant rate during the whole simulation time and an ideal-age tracer. The ideal-age tracer is defined as a tracer whose mixing ratio ξ_ia obeys the continuity equation (Waugh and Hall, 2002)

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial {\mathit{\xi }}_{\mathrm{ia}}}{\partial t}}+\mathcal{L}\left({\mathit{\xi }}_{\mathrm{ia}}\right)=\mathrm{1}\end{array}$$

(where ℒ is the advection–diffusion operator), and boundary condition ξ_ia=0 at the surface. The ideal-age tracer is transported as a regular gaseous tracer and updated at every model time step Δt with the unity tracer correction:

$$\begin{array}{}\text{(13)}& M_{\mathrm{ia}}\mapsto \left\{\begin{array}{ll}\mathrm{0},& \text{at lowest layer},\\ M_{\mathrm{ia}}+M_{\mathrm{unity}}\mathrm{\Delta }t& \mathrm{otherwise},\end{array}\right.\end{array}$$

where M_ia and M_unity are masses of the ideal-age tracer and of the unity tracer in the grid cell. The mixing ratio of the ideal-age tracer is a direct measure of the mean age of air in a cell, so the tracer is a direct Eulerian analogue of the time-tagged Lagrangian particles with clock reset at the surface. Note that the AoA derived from the ideal-age tracer and AoA from a passive tracer with a linearly growing near-surface mixing ratio are equivalent (Waugh and Hall, 2002), and implementation of both provides a redundancy needed to ensure self-consistency of our results.

The simulations were performed with four eddy-diffusivity profiles described in Sect. 3.1 and the corresponding destruction rates of sf6 and sf6nograv tracers in the uppermost model layer. All runs were initialized with the mixing ratios from the final state of a special initialization run. The initialization simulation with 0.1-Kz eddy diffusivity was started from 1970 with zero fields for all tracers, except for the unity tracer that was set to unity mixing ratio. The simulation used 1970–1989 emissions for SF₆ species from the same inventory as for the main runs (Rigby et al., 2010), and it was driven with the twice repeated ERA-Interim meteorological fields for 1980–1989. The mixing ratios of all SF₆ tracers at the end of the initialization run were scaled to match the total SF₆ burden of 20.17 Gg in 1980 (Levin et al., 2010).

4 Sensitivity and validation of SF₆ simulations

4.1 Gravitational separation and mesospheric depletion

To evaluate the relative importance of gravitational separation, mesospheric depletion, and their effect on the SF₆ concentrations, we compared the simulations for the SF₆ tracers and evaluated the relative reduction of the SF₆ content in the stratosphere due to these processes. As a conservative estimate of the reduction, we evaluated the relative differences between the tracers in the latitude belt of 70–85^∘ S, since both processes have the most pronounced effect in the southern polar vortex, where the downwelling of Brewer–Dobson circulation is the strongest.

Figure 4 The relative reduction of the SF₆ content (in %) at 70–85^∘ S due to gravitational separation with (a, b) and without (c, d) depletion and due to the combined effect of depletion and separation (e, f) at two extreme Kz cases. Note the different colour scales for (e) and (f).

Hereafter we quantify the relative difference between atmospheric contents of two SF₆ tracers, “X” and “Y” as

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }(\text{“X”},\text{“Y”})=\mathrm{2}{\displaystyle \frac{{\mathit{\xi }}_{\mathrm{X}}-{\mathit{\xi }}_{\mathrm{Y}}}{{\mathit{\xi }}_{\mathrm{X}}+{\mathit{\xi }}_{\mathrm{Y}}}}\cdot \mathrm{100}\mathit{\%}.\end{array}$$

The relative differences for the SF₆ tracers in the southern polar region (70–85^∘ S) simulated with two extreme K_z profiles is given in Fig. 4 as a function of time and altitude. Note that every 5 % of the decrease of SF₆ with respect to its passive counterpart corresponds to about 1 year of a positive bias in AoA derived from the SF₆ mixing ratios.

The reduction of the SF₆ content due to gravitational separation, if the mesospheric depletion is disabled, is given by the relative difference of sf6nochem and sf6pass (Fig. 4a, b). Expectedly, the effect of gravitational separation is most pronounced for the case of low eddy diffusivity (0.001-Kz), and the reduction of SF₆ in the altitude range of 30–50 km reaches 2 %–5 %. In the case of strong mixing, the effect of separation is about 1 %.

The reduction of the SF₆ content due to gravitational separation in the presence of stratospheric depletion is given by the relative difference of sf6nograv and sf6 tracers. The effect of the separation for low K_z is very similar between the depletion and no-depletion cases (Fig. 4c vs. Fig. 4a). Depletion reduces the effect of the gravitational separation for high K_z (Fig. 4b vs. Fig. 4d). Regardless of depletion, stronger K_z reduces the effect of the gravitational separation; however, the latter is still non-negligible if precisions of the order of a month for AoA are required.

The combined effect of depletion and gravitational separation is seen in the relative difference of sf6pass and sf6 tracers (Fig. 4e and f). For both K_z cases, the effect of depletion is stronger than the diffusive separation by more than 1 order of magnitude. Regardless of the K_z profiles, the reduction exceeds 50 %, which roughly corresponds to 10 years of an offset in the apparent AoA.

In all cases the reduction of the SF₆ content has a strong annual cycle associated with the cycle of the downwelling in winter and the upwelling in summer. Besides, the reduction has a noticeable inter-annual variability that poses substantial difficulties for applying a consistent correction to the apparent AoA. Contrary to the former two comparisons, strong eddy mixing leads to a strong reduction of SF₆ since it intensifies the transport to the depletion layers and thus enhances the depletion rate.

The simulations for different K_z have been initialized with the same state obtained from a separate spin-up simulation with 0.01-Kz, which was scaled to match total burden of SF₆ in 1980. Thus a relaxation of the SF₆ vertical distribution during the first few years of the simulations is clearly seen in Fig. 4. For the 1-Kz case (Fig. 4f), the gradual increase of the difference between SF₆ and its passive version in the troposphere can be seen as well. The rate of this increase is about 0.5 % per 39 years of the simulations. This rate should not be confused with the depletion rate of SF₆ in the atmosphere since the difference is a combined effect of depletion and growth of emission rate, despite the fact that the latter is exactly the same for both tracers.

The above comparison indicates that depletion has the stronger effect on the SF₆ mixing ratio in the upper stratosphere than gravitational separation and molecular diffusion. However, the important role of molecular diffusion in the model is that it maintains the upward flux towards the mesosphere in the simulations even if the eddy diffusivity ceases.

Further in this paper only the sf6pass and sf6 tracers will be used.

4.2 Evaluation against balloon profiles

The tropospheric concentrations of SF₆ in our simulations have been insensitive to the SF₆ destruction or to the eddy-diffusivity profiles in the stratosphere. The difference in the modelled profiles can, however, be seen above the tropopause. For comparison, we took the simulations with prescribed eddy diffusivity in the stratosphere (1-Kz, 0.03-Kz, and 0.001-Kz; see Sect. 3.1) and with dynamic eddy diffusivity ECMWF-Kz. The simulations were matched with the stratospheric balloon observations (Fig. 5) published by Patra et al. (1997), Engel et al. (2006), Ray et al. (2014), and Ray et al. (2017).

Figure 5 Observed SF₆ balloon profiles and corresponding daily-mean SILAM profiles for the date of observations. The observational data obtained from Patra et al. (1997), Ray et al. (2014, 2017), and Engel et al. (2006) are for (a)–(d) correspondingly. The observation uncertainties are about 2 % (1σ) for Hyderabad profiles (a) and smaller than the size of the symbol for Kiruna profiles (b, d). The model profiles from the WACCM model are from Ray et al. (2017).

Two balloon profiles observed at Hyderabad (17.5^∘ N, 78.6^∘ E) in 1987 and 1994 by Patra et al. (1997) indicate an increase of the SF₆ content during the time between the soundings (Fig. 5a). Both profiles have a clear transition layer from tropopause at ∼17 km to the undisturbed upper stratosphere above ∼25 km. The simulated profiles agree quite well with the observed profiles, except for the most diffusive case that gave notably smoother profiles and somewhat overstated SF₆ mixing ratios due to too strong upward transport by diffusion through the tropopause and in the lower stratosphere.

The profile in Fig. 5b has been obtained from Kiruna (68^∘ N, 21^∘ E) in early spring 2000 during the SAGE III Ozone Loss and Validation Experiment, SOLVE, (Ray et al., 2002) with the lightweight airborne chromatograph (Moore et al., 2003). The profile is affected by the polar vortex and clearly indicates a strong reduction of SF₆ with height with a pronounced local minimum at 32 km. The corresponding SILAM profiles tend to overestimate the SF₆ volume mixing ratio (vmr). The SF₆ profiles for ECMWF-Kz and 0.001-Kz match each other, since vertical mixing is negligible in both cases. The most diffusive profile, 1-Kz, has the strongest depletion in the upper part but the largest deviation from the observations below 20 km. The intermediate-diffusion profile (0.03-Kz) is almost as close to the observations as the non-diffusive profile. Moreover, the 0.03-Kz profile has a minimum at the same altitude as the observed one, albeit the modelled minimum is substantially less deep.

For comparison, Fig. 5b also contains monthly-mean profiles from the WACCM simulations by Ray et al. (2017). The WACCM profiles match very well with the observations below 17 km but turn nearly constant above, thus under-representing the depletion of SF₆ inside the polar vortex. Monthly-mean SILAM profiles (not shown) were much closer to the plotted daily profiles than to the ones of WACCM. However, the WACCM simulations did not include the electron attachment mechanism.

For the mid-latitude profile in Fig. 5c from Aire-sur-l'Adour, France (43.7^∘ N, 0.3^∘ W), all SILAM profiles except for 1-Kz fall within the observational error bars provided together with the data by Ray et al. (2017). Similar to the Kiruna case in Fig. 5b, the SILAM profiles are smoother than the observed ones and are unable to reproduce the sharp transition at 20 km.

Another profile from within the polar vortex (Fig. 5d) was observed at the same Kiruna site as the one in Fig. 5b, but three years later. The observed profile also has a minimum that is much deeper than in the modelled profiles. Similar to the case in Fig. 5b, the 0.03-Kz profile is the only one that has a pronounced minimum at the same altitude as the observed one. The minimum is a result of the spring breakdown of the polar vortex when a regular downdraught ceases and atmospheric layers decouple from each other. The reduced depth of the modelled minimum is probably caused by insufficient decoupling of the layers in the driving meteorology.

In all above cases, the 1-Kz profile is clearly far too diffusive in the non-polar cases, whereas for the Kiruna cases it overstates the lower part of the profiles and smears out the vertical structure of the profiles further above the tropopause. The SF₆ profiles simulated with ECMWF-Kz and 0.001-Kz match each other in all simulations, since vertical mixing is negligible in both cases. The SF₆ resulting from the 0.03-Kz case appears to be the most realistic out of the four considered simulations: they are close to the observed ones and have the local minima at the correct altitudes for both Kiruna profiles.

4.3 Evaluation of SF₆ against MIPAS data

The MIPAS observations provide the richest observational dataset for the stratospheric SF₆ profiles. However, each individual observation has a substantial retrieval noise error, which is noticeably larger than the difference between the observation and any of the SILAM simulations. The largest diversity of the modelled SF₆ profiles was observed in polar regions; therefore, below we show the mean profiles for each season in the southern and the northern polar areas. Besides that, we consider statistics of the model performance against MIPAS measurements in the lower and upper stratosphere separately. For simplicity, we do not show the statistics for the ECMWF-Kz runs, since they are very similar to the ones for 0.001-Kz.

Figure 6 Seasonal mean co-located SILAM SF₆ and MIPAS profiles for 2007, for southern and northern polar regions. Typical ranges covering 75 % of the averaging kernel are given with the error bars at the right-hand side of each panel. The horizontal error bars indicate systematic uncertainties of the observations that are fully correlated among profiles and do not cancel out when averaging over a large number of measurements. Dashed lines are zonal-mean SILAM profiles for a given season taken without co-location.

For the comparison, the daily-mean model profiles were co-located to the observed ones in space and time, after which an averaging kernel of the corresponding MIPAS profile was applied to the SILAM profile. For the comparison, we took only the data points with all of the following criteria met:

• MIPAS visibility flag equals 1;

• MIPAS averaging kernel diagonal elements exceed 0.03;

• MIPAS retrieval vertical resolution, i.e. the full width at the half maximum of the row of the averaging kernel, is better than 20 km;

• MIPAS volume mixing ratio noise error of SF₆ is less than 3 pmol mol⁻¹.

The mean seasonal profiles of the SF₆ mixing ratio for southern and northern polar regions derived from the MIPAS observations and the SILAM simulations for 2007 are given in Fig. 6. In order to facilitate the comparison of our evaluation with the earlier study of Kovács et al. (2017), we have chosen the same year and same layout of the panels as Fig. 3 there. The main differences between Kovács et al. (2017) and the current evaluation are the following.

• We used averages of co-located model profiles (bold lines). The non-co-located seasonal- and area-mean model profiles are given as thin dashed lines for comparison.

• We use a newer version of the MIPAS SF₆ data with considerably larger values (up to 0.6 pmol mol⁻¹) in the upper stratosphere, compared to the version that was used by Kovács et al. (2017).

• The horizontal error bars for the observed data indicate that the systematic error component is fully correlated among the profiles and does not cancel out by averaging or, in other words, the estimate of a possible bias, as analysed by Stiller et al. (2008). These errors are of the order of 4 % (below 30 km) up to 10 % (at 60 km). The contribution of the retrieval noise error is essentially negligible due to averaging. The error bars shown by Kovács et al. (2017) are noticeably larger, probably indicating that they are for the individual observed values rather than the uncertainties of the mean.

• We use 3 km vertical bins for the profiles to make the points in the MIPAS profiles distinguishable.

• We also plot the vertical extent of the averaging kernels corresponding to their half widths.

First of all, there is a substantial difference between the co-located and non-co-located model profiles. The difference is caused by the uneven sampling of the atmosphere by the satellite both in space and in time. In particular, MIPAS, being a polar-orbiting instrument, makes more profiles per unit area closer to the pole than further away. The difference gets somewhat reduced if one uses equal weights for all model grid cells instead of area-weighted averaging, especially for wide latitude belts. The major difference comes probably from the inability of MIPAS to retrieve SF₆ profiles in the presence of polar stratospheric clouds that clutter lower layers of the stratosphere and make the sampling of polar regions quite uneven both in time and in the vertical. This hypothesis agrees with the fact that the difference is most pronounced for the winter pole, especially for the South Pole in JJA, and almost invisible at a summer pole.

The comparison in Fig. 6 shows that the profiles from the SILAM simulations agree quite well to the observations in the altitude range below 20–25 km, with the most diffusive, 1-Kz, slightly overestimating the SF₆ mixing ratios. In the range above 25 km, the 1-Kz profiles indicate a decrease of SF₆ with altitude that is too fast. The 0.03-Kz profiles give the best results up to ∼40 km, except for the South Pole in JJA and the North Pole in DJF.

An interesting feature of the winter-pole MIPAS profiles is an increase of the SF₆ mixing ratio above 40 km. This increase might be caused by issues with retrievals as the systematic errors of the retrievals increase with altitude. However, non-monotonic profiles can occur due to the mean atmospheric dynamics (see the non-co-located 0.001-Kz profile in Fig. 6g).

None of the model setups are capable of reproducing the observations above 40 km. Wintertime poles also pose a problem to the model. The disagreement indicates a deficiency in the model representation of air flows in the upper part of the domain caused by insufficient vertical resolution of ERA-Interim in the upper stratosphere and lower mesosphere and a lack of pole-to-pole circulation. This discrepancy is in line with the comparisons in Fig. 5 for polar regions. The model tends to overstate the SF₆ content in the lower part of the polar vortex and understate it above 40 km.

Figure 7 Time series of monthly scores for the SILAM SF₆ mixing ratios for the whole period of the MIPAS observations in the altitude range of 10–35 km. The statistics are standard deviations of model–measurement difference (SD), absolute bias, and normalized mean bias (NMB). The statistics of the model mixing ratios extracted at nominal MIPAS altitudes are given as thin lines.

Figure 8 Same as in Fig. 7 but for the MIPAS altitude range of 30–60 km.

We also computed statistical scores of the simulated SF₆ mixing ratios for each month of the MIPAS mission. The statistics were computed separately for the altitude ranges of 10–35 km (Fig. 7) and 30–60 km (Fig. 8). As the difference in the statistical scores between the three selected simulations is quite minor, we used only observations with the retrieval target noise error below 1 pmol mol⁻¹.

The root-mean-square error turned out to be mostly controlled by the bias, and it does not allow for a clear distinction between the simulated cases. In order to disentangle the effect of bias, we have calculated the standard deviation of the model–measurement difference (SD), absolute bias, and normalized mean bias (NMB):

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathrm{SD}\left(\mathrm{pmol}{\mathrm{mol}}^{-\mathrm{1}}\right)& {\displaystyle }={〈{\left(M-\langle M\rangle -O+\langle O\rangle \right)}^{\mathrm{2}}〉}^{\mathrm{1}/\mathrm{2}},\text{(16)}& {\displaystyle }\mathrm{Bias}\left(\mathrm{pmol}{\mathrm{mol}}^{-\mathrm{1}}\right)& {\displaystyle }=〈M-O〉,\text{(17)}& {\displaystyle }\mathrm{NMB}\left(\mathit{\%}\right)& {\displaystyle }=\mathrm{2}〈{\displaystyle \frac{M-O}{M+O}}〉\cdot \mathrm{100}\mathit{\%},\end{array}$$

where M and O are modelled and observed values, respectively, and 〈⋅〉 denotes averaging over the selected model–observation pairs for the given range of times and altitudes. Along with the SD, we have plotted the RMSE of the observations due to the retrieval noise in the original MIPAS data, labelled as “MIPAS noise” in the top panels of Figs. 7 and 8.

In the altitude range of 10–35 km, the SD of model–measurement difference is uniform in time with minor peaks in August–September (Fig. 7). The level of the noise error constitutes about 85 % of the total model–measurement difference. Application of the averaging kernel to the model profiles reduces the SD. The intermediate-diffusivity case, 0.03-Kz, clearly shows the least SD uniformly over the whole observation period; the same case indicates the least absolute bias.

In the range of 30–60 km altitudes (Fig. 8), the level of the retrieval noise is noticeably higher than in the lower stratosphere. The least biased case is 1-Kz, which, however, has the largest SD. The SDs of 0.03-Kz and 0.001-Kz are on par, but the latter has the strongest bias. Thus for this altitude range the intermediate-diffusivity case also shows the best performance.

Note the slight increase of the model bias after 2009, which is likely caused by our overestimating of the emission rates since that time (see Sect. 3.4). This increase of the bias does not appear in Fig. 8 due to the delay in the response of the content in the upper layers to the changes in surface emissions.

4.4 Lifetime of SF₆ in the atmosphere

In order to estimate the atmospheric lifetime of SF₆, we turned off the emission of all SF₆ tracers in July 2016 and let the model run until the end of 2018 without emissions (Fig. 9). The decrease of the simulated burden after the emission stop can be used to estimate the removal rate from the atmosphere.

Figure 9 The time series of mean mixing ratio of SF₆ in the atmosphere simulated with emissions stopped in July 2016. The total burdens by Levin et al. (2010) and by Rigby et al. (2010) are shown for comparison.

Table 1 SF₆ destruction rate after stopping the emissions and corresponding lifetimes. Mid-2011 burden of 1.27×10⁹ moles is used as a reference for the lifetime estimate.

Time series of the total burden of SF₆ in the atmosphere in the simulations are given in Fig. 9. For easier comparison to the observed mixing ratios, the burden has been normalized with 1.78×10²⁰ moles – the total amount of air in the atmosphere – to get the mean mixing ratio. The tabulated values for the atmospheric burden of SF₆ from Levin et al. (2010) and Rigby et al. (2010) are given for comparison. Since the removal of SF₆ from the atmosphere is mostly controlled by the transport towards the depletion layer, the vertical exchange is the key controlling factor.

The decrease of the atmospheric SF₆ content after the emission stop is given in the inset in Fig. 9. As expected, after July 2016 the content of passive SF₆ stays constant, while the others begin to decrease at a rate that depends on the transport properties in the stratosphere with the faster removal for the stronger eddy diffusivity. The removal rate is driven by the SF₆ content in the upper stratosphere, which is not in equilibrium with the total atmospheric content. A typical delay between the SF₆ mixing ratio in the troposphere and the upper stratosphere, i.e. the AoA in the topmost model layer, is about 5–6 years. Therefore, for a reference we used the total amount of atmospheric SF₆ 5 years before the emission stop, i.e. 1.23×10⁹ mol, which corresponds to the mean mixing ratio of 7 pmol mol⁻¹. Dividing the destruction rate with the reference amount, one gets the range of corresponding simulated SF₆ lifetimes in the atmosphere: 600 to 2900 years. Despite the range of the tested diffusivities of 3 orders of magnitude, the loss rate varies only within a factor of 5 (Table 1).

The term “lifetime” implies a linear decay; however, due to emissions the distribution of SF₆ in the atmosphere is far from equilibrium, so the decay is not proportional to the burden. A more accurate way to estimate the lifetime would be to perform a multi-decade simulation without sources to get the distribution into a quasi-equilibrium with the mesospheric sink. In such a quasi-equilibrium the model of linear decay of SF₆ in the whole atmosphere becomes applicable and the lifetime can be estimated as a simple ratio of the burden to the loss rate. The uncertainty in the equilibrium burden corresponding to the modelled loss rates in Table 1 can be estimated as the range of AoA in the upper stratosphere (∼0.5 years) divided by the growth rate of the burden (0.04 yr⁻¹), i.e. about 2 %. A larger uncertainty comes from the over-simplistic parametrization of the loss in the model, which is more difficult to quantify.

The best-performing simulation, 0.03-Kz, resulted in 1540 years lifetime. Given the uncertainties above, it meets the ranges suggested by earlier studies. It is in a good agreement with the range of 800–3200 years from the model studies (Ravishankara et al., 1993; Morris et al., 1995), and it is close to the upper bound of the 580–1400 years range recently obtained by Ray et al. (2017) from the balloon profile given in Fig. 5b.

Our estimate is also slightly above the range given by Kovács et al. (2017), who obtained 1120–1475 years. However, in the simulations by Kovács et al. (2017) the mixing ratios of SF₆ in the stratosphere and the lower mesosphere were noticeably higher than those retrieved by MIPAS and practically flat in the range of 30–50 km. Such modelled profiles likely indicate a vertical exchange in the model that is too strong; a loss that is too strong, as a consequence; and corresponding low bias of the estimated lifetime.

5 Simulations of AoA

5.1 Eddy diffusivity and simulated AoA

The effect of the vertical eddy diffusivity on AoA in the stratosphere was evaluated with the same set of three prescribed and one dynamic K_z profiles, as for SF₆ simulations. An example of annual-mean distributions of AoA is given in Fig. 10. The Hunten (1975) K_z profile (Fig. 10a) gives AoA in the stratosphere of about 3.5 years. It is much shorter than the estimates of the stratospheric AoA (e.g. Waugh, 2009; Engel et al., 2009) from the observations of various tracers. Three other profiles of K_z result in almost identical average distributions of AoA with typical stratospheric AoA of 5.5 years, which agrees quite well with the experimental estimates. In these cases AoA is controlled by the transport with mean winds. Since 0.03-Kz profiles result in the most realistic distribution of SF₆ in our simulations, in the current section we will use simulated distributions of tracers with this parameterization.

Figure 10 The zonal-mean spatial distribution of the ideal-age AoA for 2011 calculated for different eddy-diffusivity profiles.

5.2 AoA and apparent SF₆ AoA

The AoA for all tracers (except for the ideal age) was calculated as a simple time lag between the mixing ratio at each point of the domain and the mean near-surface mixing ratio. As it has been pointed out by Waugh and Hall (2002), this lag equals to AoA only in the case of a fully passive tracer with linearly growing (or decreasing) near-surface mixing ratio. Corrections have been applied to the AoA derived from SF₆ in many studies (Volk et al., 1997; Stiller et al., 2008, 2012; Engel et al., 2009) to account for non-linear growth of the near-surface SF₆ mixing ratio and for its mesospheric sink. The corrections rely heavily on various assumptions that can hardly be rigorously verified. Therefore, in this study we do not apply any corrections to the AoA derived from the time lags of tracers. The corrections and assumptions behind them are discussed in Sect. 6.

Figure 11 Zonal-mean distributions of atmospheric AoA simulated with passive, ideal-age, and two SF₆ tracers (average for 2012).

The constant-rate emission of the passive tracer resulted in almost linear growth of its near-surface mixing ratio after the spin-up. The latter makes the age derived from the passive tracer equivalent to the age derived from the ideal-age tracer. The resulting distributions are indeed very close to each other (Fig. 11a and b). The agreement confirms the self-consistency of the transport procedure since the tracers have opposite sensitivity to the advection errors: higher mixing ratios correspond to younger air for the accumulating tracers, while for the ideal-age tracer higher mixing ratios correspond to older air. The remaining differences are caused by spatial inhomogeneities of near-surface mixing ratio of the passive tracer due to variations in the near-surface air density.

The distribution of the AoA derived from sf6pass (Fig. 11c) is qualitatively similar to the ideal-age one; however, one can see substantial differences. The negative AoA in the northern troposphere for the sf6pass tracer is caused by the predominant location of the sources in the Northern Hemisphere, so the concentrations there exceed the global-mean levels. The growing rate of the SF₆ emissions leads to the faster-than-linear increase of near-surface mixing ratios, which leads to an old bias of up to 3–5 months of the sf6pass AoA. This old bias has been one of the drawbacks of the SF₆ AoA pointed out by Garcia et al. (2011).

The ages shown in Fig. 11a–c agree well with the ages derived from the in situ observations of SF₆ and CO₂ at the 25 km altitude by Waugh and Hall (2002). They also agree quite well with the earlier simulations with five climate models that give annual mean ages in the upper stratosphere between 4.5 and 5.5 years (Butchart et al., 2010). The simulations result in about 1–1.5 years younger air than diabatic mean age obtained with the Lagrangian model computations of Diallo et al. (2012) (Fig. 11 is directly comparable with Fig. 2 there) and about 1 year older air than kinematic mean age. Since our preprocessor of wind fields differed strongly from that by Diallo et al. (2012), this similarity is an important indicator of consistency of the numerical procedures applied in both studies.

A substantial disagreement, however, exists with the ages derived from the MIPAS satellite observations (Stiller et al., 2012; Haenel et al., 2015). The authors calculated ages exceeding 10 years in the polar areas and in the upper stratosphere. The reason for the disagreement follows from the above analysis: SF₆ can neither be considered a passive tracer nor does its mixing ratio in the troposphere grow linearly with time. Denoting the AoA derived from the SF₆ profiles as “apparent AoA” (Waugh and Hall, 2002), we calculated it from the SILAM-predicted SF₆ profiles, which, as shown above, agree well with AoA derived from MIPAS. The resulting model-based apparent AoA (Fig. 11d) is indeed much older than the ideal-age AoA. The distribution of the apparent SF₆ AoA agrees with the AoA retrieved from MIPAS SF₆ profiles by Haenel et al. (2015): well over 5 years AoA around the Equator with well over 10 years AoA in the polar regions.

The effect of the apparent over-ageing in the stratosphere due to the subsidence of the mesospheric air was estimated by Stiller et al. (2012) to be a fraction of a year in the upper stratosphere. Earlier experimental balloon studies (Strunk et al., 2000) indicated an up to 3.5-year difference between CO₂ and SF₆ ages. In our simulations, the over-ageing due to the SF₆ depletion and other factors discussed in the previous sections is much stronger and affects the whole stratosphere.

5.3 Trends in apparent AoA

Changes in the AoA have been used in many studies as an indicator of changes in the atmospheric circulation. In order to evaluate the effect of the way the AoA is computed on its trend, we have calculated trends of the apparent AoA at different altitudes and latitudes for 11 years (2002–2012). This period roughly covers the MIPAS mission and allows for comparison with trends reported by Haenel et al. (2015).

Figure 12 Vertical profiles of the simulated age of air linear trends over 11 years (2002–2012) for example latitude belts.

Figure 13 Anomaly of the ideal-age AoA (years) for the period of 1990–2018 with respect to the mean AoA.

The zonal-mean vertical profiles of the AoA trends during 2002–2012 are shown in Fig. 12 for five latitudinal belts. The presented variable is a slope of the linear fit of the deseasonalized monthly-mean time series for each tracer, averaged over the corresponding latitudinal belt and the model layer. The fit was made with the ordinary least-squares method. The error bars show 95 % confidence intervals calculated as if a model of linear trend with uncorrelated Gaussian noise was applicable to the time series.

The trends of the apparent AoA for the non-passive SF₆ species have a clear increase with height in the upper part of the profiles. Such behaviour agrees well with the AoA trends by Haenel et al. (2015, Fig. 7) obtained from the MIPAS observations. The over-ageing due to the mesospheric depletion of SF₆ has been discussed and estimated by Haenel et al. (2015) and Kovács et al. (2017). However, Fig. 12 shows that the mesospheric depletion of SF₆ also affects its trend: the over-ageing increases with time. The reason is that depletion is proportional to the SF₆ load, which grows with time. This effect has been pointed out and evaluated earlier for N₂O by Schoeberl et al. (2000). For SF₆, the effect of its loss on the AoA was evaluated by Stiller et al. (2012), who concluded that “in-mixing of mesospheric SF₆-depleted air plays a minor role for the assessment of AoA trends”, at least within the framework of their approach (2002–2010, up to 35 km altitude).

The apparent AoA derived with the passive SF₆ tracer sf6pass indicates a negative trend of about 0.5 years per decade. The trend is caused by the temporal variation of SF₆ emissions. In order to get an unbiased AoA estimate from the passive tracer, one needs the mixing ratio at the surface to be increasing linearly with time. A steady growth of emission rate leads to the faster-than-linear increase of the near-surface mixing ratio and thus a low bias of the AoA. According to the inventory (Levin et al., 2010) used in this study, the SF₆ emission rate was growing in 1997–2000 about twice slower than after 2005. Consequently, the negative bias of the apparent AoA has increased resulting in the negative trend of the AoA in the stratosphere.

The AoA trends derived from the ideal-age and passive tracers agree through the whole range of altitudes and latitudes, indicating internal consistency of our simulations. The main common feature of the profiles is the negative tendency of about −0.5 year per decade in the altitude range of 15–30 km with a profile that varies across altitudes. Similar-magnitude trends for the same period were reported by Plöger et al. (2015), who used the same ERA-Interim to simulate AoA. The major difference between the obtained trends is that we have consistently negative trends for both hemispheres, whereas Plöger et al. (2015) indicate a positive trend as a fraction of year per decade in the altitude range of 20–30 km in the Northern Hemisphere and a similar-magnitude negative trend in the Southern Hemisphere. The reason for the discrepancy deserves further investigation. A possible reason for the discrepancy is that Plöger et al. (2015) used diabatic heating rates as vertical velocity, and it is known that the diabatic and kinematic vertical transport is inconsistent in the reanalysis (Abalos et al., 2015).

The trends might be a feature of the non-uniformity of the ERA-Interim dataset, which was produced with assimilation of an inhomogeneous set of the observations. During 2002–2012, the amount of the assimilated data of the upper-air temperatures was an order of magnitude higher than before 2000 and 2 orders of magnitude higher than after 2010 (Dee et al., 2011). It had a clear impact on the patterns of the analysis increments in ERA-Interim and, consequently, on the predicted stratospheric circulation. Due to such inhomogeneities, the quality of trends derived from reanalysis data needs to be verified for each geophysical quantity (Dee et al., 2011). Deducing reliable trends for atmospheric temperature, a quantity that is measurable and extensively assimilated, took a major effort (Simmons et al., 2014). The fact that the AoA is not a directly observable quantity makes the verification of the AoA trends difficult.

To get more insight into the nature of the simulated long-term AoA variability at different altitudes and latitudes, we have plotted the time series of the monthly zonal-mean ideal-age AoA for the same latitude belts as in Fig. 12 over 1990–2018 (Fig. 13). To make the temporal variations more visible, the mean AoA profile for each latitude averaged over the same period was subtracted from the profiles. One can see a clear seasonal variation of the AoA outside the equatorial zone. The variation has opposite phases in the upper and the lower stratosphere. In the altitude range of 20–30 km, where the trends are most pronounced, the temporal variation of the AoA has a ramp structure with more-or-less steady intervals and relatively quick changes. Such a structure is similar to the one shown for the ERA-Interim analysis increments (Dee et al., 2011) and is likely to be caused by temporal inhomogeneities in the assimilated dataset. Therefore, we do not draw any conclusion here on the actual trends of AoA, but we highlight that trends of the apparent AoA are strongly influenced by the selected time interval and by the method of the trends calculation.

6 Discussion

The present study has several limitations that deserve specific attention. Forced zero air flux through the domain top at 0.1 hPa caused distortion of the mean transport within the domain and left diffusive transport as the only means for the upper-boundary fluxes of SF₆. Secondly, we used prescribed profiles of the eddy diffusivity within the domain, which also affects the results of the simulations. In this section we evaluate the role of these distortions.

6.1 Distortions of air flows

The transport procedure used in this study is done with a “hardtop” diagnostics, forcing zero mass fluxes at the domain top and forced air-mass conservation everywhere within the domain. Since the upper boundary of the domain is at 0.1 hPa, the divergence of the air flow above that level in the meteorological data used to drive the model is compensated by adjusting the divergences within the domain. To evaluate the effect of this adjustment on the mean circulations, we used the new ERA-5 dataset, which has the topmost level at 10⁻³ hPa. The diagnostic procedure was applied to ERA5 for two sets of vertical layers: the 61 ERA-Interim layers, same as used in the SF₆ simulations (hereafter ERA5-cut), and a refined vertical matching the 137 native ERA5 vertical layers (hereafter ERA5). The resulting vertical winds were compared to the ones used in the SF₆ simulations: 61 layers diagnosed from ERA-Interim. The seasonal and zonal-mean vertical air-mass fluxes, expressed in units of pascals per day (Pa d⁻¹), for the three cases and two solstice seasons of 2017, are shown in Fig. 14 together with the corresponding layer boundaries.

Figure 14 The seasonal and zonal-mean vertical air-mass fluxes diagnosed by SILAM from ERA5 and ERA-Interim fields for 2017 solstice seasons, expressed in terms of vertical velocity ω. Updraughts are red. The vertical-layer boundaries are shown with grey lines.

The wind patterns in ERA5 (Fig. 14a, b, d, e) have finer features than in ERA-Interim due to the higher horizontal resolution. The difference between the ERA5 and ERA5-cut vertical winds is the strongest at the cut-domain top (0.1 hPa, 65 km), where the zero vertical air-mass flux is forced. For both seasons, the disturbances introduced by the cut vertical are minor, except for the summertime poles (South Pole in Fig. 14a, b and North Pole in Fig. 14d, e), where a noticeable disturbance is visible down to 35–40 km altitude. Such systematic disturbances influence the performance of the AoA and the SF₆ simulations in the polar stratosphere, and they are a probable reason for the failure of the model to reproduce the SF₆ profiles there (see Fig. 6).

The comparison of the mass fluxes for the same vertical levels (panels b vs. c or e vs. f in Fig. 14) shows that the difference between ERA-Interim and ERA5 is noticeably larger than between cut and full vertical of ERA5. Thus we conclude that the distortions introduced by our diagnostic procedure are within the uncertainty of the input meteorological data.

6.2 Top-boundary mass fluxes and eddy diffusion profiles

The used modelling approach replaces the vertical transport through the domain top with the diffusive fluxes for the depleting SF₆ and a hard lid for other species. This approach is unlikely to introduce major disturbances into the AoA fields since the AoA is quite uniform close to the domain top. The uncertainty introduced with this approach into the SF₆ fields is not straightforward to evaluate due to a major uncertainty in the vertical diffusivity profiles.

As mentioned in Sect. 3.1, the eddy-diffusivity profiles of the C-IFS model from the ERA5 reanalysis (Fig. 1) are clearly unrealistic within and above the stratosphere. They do not exhibit any growth of the eddy diffusivity in the mesosphere either. According to Lindzen (1981) the mean diffusivity due to the breaking gravity waves has an order of magnitude of 10² m² s⁻¹, whereas the eddy diffusion in ERA5 for that region is below the molecular diffusivity (Fig. 1). On the other hand, if we assume that the mesospheric turbulence results in a diffusivity profile as predicted by Lindzen (1981) (Fig. 1), then such turbulence provides quite rapid exchange of SF₆ towards the depletion layers, making the advective vertical transport above ∼50 km negligible. The profiles of Lindzen (1981), however, do not allow for a simple extrapolation below 50 km; therefore, the vertical profiles by Massie and Hunten (1981) (1-Kz) were involved as the ones that are simple to implement and smooth enough to be easily approximated and extrapolated.

The normalized diffusive SF₆ mass fluxes above the domain top for the scaled profiles of the eddy diffusivity (Fig. 3) allow for evaluation of the validity of the assumption of neglected regular vertical transport above the domain top. The equivalent vertical air-mass flux due to diffusion at the level of 0.1 hPa (domain top) is $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{6}}$, $\mathrm{9}\times {\mathrm{10}}^{-\mathrm{7}}$, and it is $\mathrm{2.5}\times {\mathrm{10}}^{-\mathrm{7}}$ kg m⁻² s⁻¹ for 1-Kz, 0.03-Kz, and 0.001-Kz. These mass fluxes, divided by g, give the vertical velocities of −5, −0.8, and −0.4 Pa d⁻¹. Comparing these values to those shown in Fig. 14 for the level of 65 km, one can see that the diffusive limit is valid for the 1-Kz profile except for the vicinities of the poles. For lower values of the eddy diffusivity, the regular circulation becomes comparable with the diffusion or even exceeds it.

Although the 0.03-Kz profiles gave better agreement with the observations of SF₆, this does not indicate that 0.03-Kz profiles are more realistic. This profile is likely to over-mix the lower stratosphere and under-mix the upper stratosphere and the mesosphere. Thus the vertical structure of the eddy diffusivity remains a major source of uncertainty in the modelling approach. Using more realistic vertical diffusion profiles and high-top ERA5 reanalysis is planned for the future studies.

6.3 Notes on the observed SF₆ age

There are three main factors responsible for the SF₆ age being different from the ideal age: the non-linear growth of tropospheric burden, the gravitational separation, and the mesospheric sink. Here we consider the effects of these factors and corrections to the SF₆ observations that can be applied to compensate for the effect of these factors on the resulting AoA.

The correction for the non-linear growth rate introduced by Volk et al. (1997) and used in many subsequent studies is based on a simple analytical model of 1D diffusion with constant diffusivity and exponential distribution of air density. The model was suggested by Hall and Plumb (1994) as an illustration for the concept of the age spectrum. The model spectrum has two parameters: the mean age Γ and the width parameter Δ. In order to use the spectrum for the correction, one has to involve an additional constraint connecting these parameters. Based on a 3D simulation with a general circulation model (GCM), Hall and Plumb (1994) suggested that a constant ratio ${\mathrm{\Delta }}^{\mathrm{2}}/\mathrm{\Gamma }=\mathrm{0.7}$ years can be used throughout the stratosphere. Note that this dimensional parameter, while having proper units originally, appears without units in several subsequent papers (Engel et al., 2002; Stiller et al., 2012). Volk et al. (1997) used the value ${\mathrm{\Delta }}^{\mathrm{2}}/\mathrm{\Gamma }=(\mathrm{1.25}\pm \mathrm{0.5})$ years for the lower stratosphere based on the results of a more advanced GCM than the one used by Hall and Plumb (1994). With this approach Volk et al. (1997) obtained the difference between the mean age and the lag time (apparent SF₆ age). The difference becomes significant for the air older than 3–4 years and approaches 0.5±0.25 years for the oldest (6 years) air, which agrees quite well with the difference between the ideal age and the passive SF₆ in our simulations (Fig. 11b, c). The correction for this difference derived from the 1D model has been used to reduce the systematic biases from the SF₆-based AoA, though “the global stratosphere is poorly represented by a 1D model” (Waugh and Hall, 2002). The uncertainty of the correction of up to ±0.5 years is systematic, is not guaranteed to be uniform in space or in time, and is likely to affect the trend estimates.

As shown in Sect. 4.1, the biases introduced to the SF₆-based AoA by gravitational separation reach a fraction of a year in the upper stratosphere. One could in principle elaborate a correction for gravitational separation; however, the correction would be well within the uncertainty of the correction for the non-linear growth rate, and thus it is probably not worth considering.

The mesospheric sink has the largest impact on the SF₆-derived AoA. The effect of the mesospheric sink is clearly visible above 15–20 km at all latitudes (Fig. 11) and leads to a strong over-ageing in the upper layers, especially in the polar areas. The effect of the sink alone can explain the discrepancy between the AoA derived from the MIPAS observations (Haenel et al., 2015) and the AoA from the modelling studies (e.g. Diallo et al., 2012; Brinkop and Jöckel, 2019). Compensating for such over-ageing is hardly possible without detailed modelling of the physical processes including depletion, diffusion, and mean transport. Since the AoA is derived as a difference of the SF₆ mixing ratios, whereas depletion introduces multiplicative change to the SF₆ abundance, the effect of the sink on apparent SF₆ AoA is unsteady in time (Fig. 12).

Once one has a model that is capable of reproducing the processes behind the SF₆ depletion, it is natural to validate such a model directly against the available SF₆ observations rather than deriving the AoA from the SF₆ observations and comparing it against the modelled one. In any case the AoA derived from the SF₆ tracer observations with all the corrections can not be considered a purely observed one.

7 Conclusions

Eulerian simulations of the tropospheric and stratospheric transport of several tracers were performed with the SILAM model driven by the ERA-Interim reanalysis for 1980–2018. The simulations included species representing SF₆ under different assumptions: a passive tracer emitted uniformly at the surface and an ideal-age tracer directly comparable to other state-of-the-art CTM simulations of the AoA. To our best knowledge this is the first systematic evaluation of AoA derived from several different tracers within the same multi-decadal simulation, combined with the extensive evaluation against MIPAS and balloon SF₆ observations.

Due to the limited vertical coverage and resolution of ERA-Interim in the upper stratosphere, the SILAM simulation domain had a lid at 0.1 hPa, which is below the altitude of the SF₆ destruction. In order to perform realistic simulations of SF₆ in our setup, the eddy diffusion in the upper stratosphere and lower mesosphere had to be parameterized, along with the mesospheric sink of SF₆.

A set of simulations with different parameterizations for the vertical eddy diffusion showed that published profiles derived with no account for advection (e.g. Massie and Hunten, 1981, and references threrin) overestimate the eddy diffusivity. On other hand, the eddy-diffusivity profiles for scalars calculated from the ERA-Interim fields, according to the IFS procedures (ECMWF, 2015) or readily available from the ERA5 reanalysis, appear to be of no relevance for the upper stratosphere, since they fall below the molecular diffusivity. Evaluation of our simulations against satellite and balloon observations indicated that the best agreement between the simulated and observed SF₆ mixing ratios within the model domain is achieved for the tabulated eddy-diffusivity profile of Hunten (1975) scaled down with a factor of 30. However, this conclusion is likely to be a feature of the specific model setup. Thus, the question of the importance and magnitude of the eddy diffusivity in the upper stratosphere and lower mesosphere remains open, and the SF₆ observations are potentially a good means of its evaluation.

The mesospheric sink of SF₆ has a major impact on the mixing ratios above 20 km. The depletion impact is especially strong in the wintertime polar areas due to the descent within a polar vortex. A set of sensitivity tests showed that molecular diffusion and gravitational separation of SF₆ are responsible for up to a few percent of further reduction in SF₆ mixing ratios in the upper stratosphere.

A good agreement of the simulated SF₆ distribution with the MIPAS observations up to the altitudes of 30–35 km and with the available balloon profiles was shown. The standard deviation between the MIPAS and the modelled SF₆ mixing ratios is mainly 80 %, controlled by the noise error of the satellite retrievals; i.e. the standard deviation between model and MIPAS is about as large as the error in the satellite data. The results of the comparison also underline the importance of accurate co-location of the modelled and the observed data in terms of space, time, and vertical averaging of the observed data.

The lifetime of SF₆ in the atmosphere estimated from the best-performing setup is about 1500 years, which is at the high side of the range of other recent estimates. Our estimate is likely to be biased high due to underrepresented vertical exchange at the domain top due to missing advective transport and the missing effect of breaking gravity waves.

Our simulations were able to reproduce both AoA obtained in other model studies and the apparent SF₆ AoA derived from the MIPAS observations. This highlights the role of fast mesospheric destruction of SF₆ due to the electron attachment mechanism. Having all tracers within the same simulations, we were able to trace the differences in the estimated AoA to the peculiarities of each tracer. A good agreement between the passive tracer AoA and the ideal-age AoA indicates a consistency of the simulations, since these two methods have opposite sign of sensitivity to errors in the transport scheme.

The mesospheric sink has severe implications for the AoA derived from the SF₆. The apparent over-ageing introduced by the sink is large and variable in space and time. Moreover, the over-ageing due to the sink increases as the atmospheric burden of SF₆ grows. All this makes SF₆ unsuitable to infer AoA above ∼20 km. For a fully passive SF₆ tracer, the variable rate of emissions causes deviations from the ideal age, and these deviations can be compensated to some extent. However, correcting the deviations due to the mesospheric sink of SF₆ is hardly possible. These deviations appear as long-term trends in the apparent AoA. These trends differ from the trends in the ideal-age AoA and have no direct correspondence to the actual trends in the atmospheric circulation.

Procedures used to derive the AoA from observations of various tracers in the atmosphere are inevitably based on assumptions and idealizations that have limited and often unknown area of applicability. The resulting uncertainties in the AoA are large enough to preclude the use of apparent AoA and its trends for evaluation of changes in atmospheric circulation or for validation of atmospheric models. Observations of the tracers themselves, however, have well quantified uncertainties, so direct comparisons of simulated tracers to the observed ones are a very promising means for the atmospheric model evaluation. AoA in turn is a convenient means for model inter-comparison if a protocol of the AoA derivation is well specified.