Three-dimensional simulation of stratospheric gravitational separation using the NIES global atmospheric tracer transport model

A three-dimensional simulation of gravitational separation, defined as the process of atmospheric molecule separation under gravity according to their molar masses, is performed for the first time in the upper troposphere and lower stratosphere. We analyze distributions of two isotopes with a small difference in molecular mass (¹³C¹⁶O₂ (M_i=45) and ¹²C¹⁶O₂ (M_i=44)) simulated by the National Institute for Environmental Studies (NIES) chemical transport model (TM) with a parameterization of molecular diffusion. The NIES model employs global reanalysis and an isentropic vertical coordinate and uses optimized CO₂ fluxes. The applicability of the NIES TM to the modeling of gravitational separation is demonstrated by a comparison with measurements recorded by high-precision cryogenic balloon-borne samplers in the lower stratosphere. We investigate the processes affecting the seasonality of gravitational separation and examine the age of air derived from the tracer distributions modeled by the NIES TM. We find a strong relationship between age of air and gravitational separation for the main climatic zones. The advantages and limitations of using age of air and gravitational separation as indicators of the variability in the stratosphere circulation are discussed.

1 Introduction

Chapman and Milne (1920) proposed two different dynamical regimes: the homosphere (perfect mixing) in the lower part of the atmosphere and the heterosphere (diffusive equilibrium), where molecular diffusion leads to a separation of the atmospheric constituents according to their molecular mass. The process of atmospheric molecule separation according to their molar mass under gravity is termed “gravitational separation” (GS). The GS of atmospheric components is recognized as dominant in the atmosphere above a level of about 100 km called the turbopause (Warneck and Williams, 2012). Recently, the existence of GS of the major atmospheric components in the stratosphere was confirmed both experimentally using a precise cryogenic sampler and theoretically by two-dimensional numerical model simulations (Ishidoya et al., 2008, 2013).

The Brewer–Dobson circulation (BDC) is the global circulation in the stratosphere, consisting of air masses that rise across the tropical tropopause and then move poleward and descend into the extratropical troposphere (Butchart, 2014). Currently, the “mean age of air” (mean AoA), defined as the average time that an air parcel has spent in the stratosphere, is perhaps best known for being a proxy for the rate of the stratospheric mean meridional circulation and the whole BDC. The mean AoA provides an integrated measure of the net effect of all transport mechanisms on the tracers and air mass fluxes between the troposphere and stratosphere (Hall and Plumb, 1994; Waugh and Hall, 2002).

The increase in greenhouse gas abundance is leading to increased radiative forcing and therefore to warming of the troposphere and cooling of the stratosphere (Field et al., 2014). A strengthened BDC under climate change in the middle and lower stratosphere is robustly predicted by various chemistry–climate models (Butchart et al., 2006; Butchart, 2014; Garcia and Randel, 2008; Li et al., 2008; Stiller et al., 2012; Garfinkel et al., 2017). However, this is not in agreement with over 30 years of observations of the age of stratospheric air (Engel et al., 2009; Hegglin et al., 2014; Ray et al., 2014).

The sparse limited observations in the upper troposphere and lower stratosphere (UTLS) mean that chemical transport models (CTMs) are complementary and useful tools for widely diagnosing the BDC and representing the global transport and distribution of long-lived species. CTMs perform relatively well in the UTLS despite resolution issues. Confidence is high in the ability of models to reproduce many of the features, including the basic dynamics of the stratospheric BDC and the tropospheric baroclinic general circulation in the extratropics, the tropopause inversion layer, the large-scale zonal mean, and tropical and extratropical tropopause (Hegglin et al., 2010; Gettelman et al., 2011).

As future changes to the BDC are likely to be complex, a suite of methods, parameters, and tools is necessary to detect these changes. Leedham Elvidge et al. (2018) evaluated the capability of using seven trace gases to estimate stratospheric mean ages. Ishidoya et al. (2013) proposed using GS as an indicator of changes in the atmospheric circulation in the stratosphere. Analyses of GS, in addition to the CO₂ and SF₆ ages, may be useful for providing information on stratospheric circulation.

Ishidoya et al. (2013) also performed the first simulation of GS using the NCAR two-dimensional (2-D) SOCRATES model (Simulation Of Chemistry, Radiation, and Transport of Environmentally important Species, Huang et al., 1998), which is an interactive chemical, dynamical, and radiative model. The spatial domain of the model extends from the surface to 120 km in altitude. The vertical and horizontal resolutions are 1 km and 5^∘, respectively. To extend the earlier work and overcome the inherent limitations of the 2-D model, we here present a more quantitative analysis using the National Institute for Environmental Studies (NIES) Eulerian three-dimensional (3-D) transport model (TM). The remainder of this paper is organized as follows. Overviews of the NIES TM, a method for modeling GS, and the simulation setup are provided in Sect. 2. In Sect. 3 we study modeled GS and compare vertical profiles with those observed. Finally, a summary and conclusions are provided in Sect. 4.

2 Model and method

For further investigation of the GS process we redesigned and modified the NIES TM, which has previously been used to study the seasonal and inter-annual variability of greenhouse gases (i.e., CO₂ and CH₄ by Belikov et al., 2013b).

2.1 Model

The NIES model is an offline transport model driven by the Japanese Meteorological Agency Climate Data Assimilation System (JCDAS) datasets (Onogi et al., 2007). It employs a reduced horizontal latitude–longitude grid with a spatial resolution of $\mathrm{2.5}{}^{\circ }\times \mathrm{2.5}{}^{\circ }$ near the equator (Maksyutov and Inoue, 2000) and a flexible hybrid sigma–isentropic (σ–θ) vertical coordinate, which includes 32 levels from the surface up to 5 hPa (Belikov et al., 2013b).

The model uses a revised version of hybrid isentropic grid. The original version used isentropic levels above the 350 K potential temperature level and sigma levels between the surface and 350 K level (Belikov et al., 2013b). A modified hybrid isentropic grid was introduced to simulate better vertical transport above the tropopause, extending the bottom level of the isentropic part to 295 K, as first used in the NIES TM simulation for the age-of-air intercomparison study (Krol et al., 2018). To limit application of the isentropic grid to the mid- to upper troposphere, for each potential temperature level between 295 and 350 K, a corresponding upper limit for pressure is set at fixed sigma level. For each theta level in a list [295, 300, 305, 310, 315, 320, 330, 340], the upper sigma limit is gradually changing from a sigma level of 0.6 for the 295 K level to 0.35 for the 340 K level, as [0.6, 0.54, 0.5, 0.47, 0.44, 0.41, 0.38, 0.35, 0.32]. For model levels between 295 and 340 K, once the sigma level reaches the prescribed maximum for this model level, vertical transport switches from one based on the diabatic heating rate to using vertical wind provided by reanalysis. Over the isentropic part of the grid, the vertical transport follows the seasonally varying climatological diabatic heating rate derived from reanalysis.

Following the approach by Hack et al. (1993), transport processes in the planetary boundary layer (height provided by the ECMWF ERA-Interim reanalysis; Dee et al., 2011) and in the free troposphere are separated with a turbulent diffusivity parameterization. The modified cumulus convection parameterization scheme computes the vertical mass fluxes in a cumulus cell using conservation of moisture derived from a distribution of convective precipitation in the reanalysis dataset (Austin and Houze Jr., 1973; Belikov et al., 2013a). To set cloud top and cloud bottom height a modified Kuo-type parameterization scheme (Grell et al., 1995) is used. Computation of entrainment and detrainment processes accompanying the transport by convective updrafts and downdrafts is as employed by Tiedtke (1989).

2.2 Molecular diffusion

According to Banks and Kockarts (1973), assuming a neutral gas, the equation for the vertical component of the diffusion velocity of gas 1 relative to gas 2 ($w=w_{\mathrm{1}}-w_{\mathrm{2}}$) in a binary mixture of gas 1 and gas 2 in a gravitational field can be written as

$$\begin{array}{}\text{(1)}& w_{\mathrm{1}}-w_{\mathrm{2}}=-D_{\mathrm{12}}[{\displaystyle \frac{n^{\mathrm{2}}}{n_{\mathrm{1}}{n}_{\mathrm{2}}}}{\displaystyle \frac{\partial (n_{\mathrm{1}}/n)}{\partial z}}+{\displaystyle \frac{M_{\mathrm{2}}-M_{\mathrm{1}}}{M}}{\displaystyle \frac{\mathrm{1}}{p}}{\displaystyle \frac{\partial p}{\partial z}}+{\displaystyle \frac{{\mathit{\alpha }}_T}{T}}{\displaystyle \frac{\partial T}{\partial z}}],\end{array}$$

where n₁ and n₂ are the concentrations of particles 1 and 2, respectively; $n=n_{\mathrm{1}}+n_{\mathrm{2}}$ is the total concentration of the binary mixture; T is the absolute neutral gas temperature; p is the total pressure; M₁ and M₂ are the masses of particles 1 and 2; $M=(n_{\mathrm{1}}{M}_{\mathrm{1}}+n_{\mathrm{2}}{M}_{\mathrm{2}})/(n_{\mathrm{1}}+n_{\mathrm{2}})$ is the mean molecular mass; D₁₂ is the molecular diffusion coefficient of gas 1 in gas 2; and α_T is the thermal diffusion factor. The three terms in the brackets represent, from left to right, the effects of concentration gradient, pressure gradient, and temperature gradient, respectively, on the molecular diffusion velocity.

If component i is a multicomponent mix of minor constituents, then, using the hydrostatic equation, the perfect gas law, and Eq. (1), the diffusion velocity w_i for the ith minor constituent can be written in time-independent form (for more details see pp. 33–34 in Banks and Kockarts, 1973):

$$\begin{array}{}\text{(2)}& w_i=-D_i[{\displaystyle \frac{\mathrm{1}}{n_i}}{\displaystyle \frac{\partial n_i}{\partial z}}+{\displaystyle \frac{\mathrm{1}}{H_i}}+(\mathrm{1}+{\mathit{\alpha }}_{Ti}\left){\displaystyle \frac{\mathrm{1}}{T}}{\displaystyle \frac{\partial T}{\partial z}}\right],\end{array}$$

where $H_i=kT/M_{i}g$ is the scale height, $k=\mathrm{1.38}\times {\mathrm{10}}^{-\mathrm{23}}$ J K⁻¹ is the Boltzmann constant, g=9.81 m s⁻² is the standard acceleration due to gravity for Earth, and D_i is the molecular diffusion coefficient. Here also n_i and ${\mathit{\alpha }}_{T_i}$ are the number density and the thermal diffusion factor for species i, respectively.

Similar to the SOCRATES model the vertical component of velocity is converted to the vertical component of the molecular diffusion flux of a minor constituent i relative to air (Huang et al., 1998):

$$\begin{array}{}\text{(3)}& f_i=-D_i[{\displaystyle \frac{\partial n_i}{\partial z}}+{\displaystyle \frac{n_i}{H_i}}+(\mathrm{1}+{\mathit{\alpha }}_{Ti}\left){\displaystyle \frac{n_i}{T}}{\displaystyle \frac{\partial T}{\partial z}}\right].\end{array}$$

The molecular diffusion coefficient is estimated from kinetic gas theory by Banks and Kockarts (1973) to be

$$\begin{array}{}\text{(4)}& D_i\left[{\mathrm{cm}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}\right]=\mathrm{1.52}\times {\mathrm{10}}^{\mathrm{18}}\sqrt{{\displaystyle \frac{\mathrm{1}}{M_i}}+{\displaystyle \frac{\mathrm{1}}{M}}}\times {\displaystyle \frac{\sqrt{T}}{n}},\end{array}$$

where M_i and M are the mass of the minor constituent i and the mean molecular mass in atomic mass units, respectively, and n is the number density of air.

To derive a diffusive flux formulation consistent with the NIES model transport equation, the number density n_i is substituted by the mixing ratio $C_i=n_i/n$:

$$\begin{array}{}\text{(5)}& f_i=-D_i\times n\times C_i[{\displaystyle \frac{\mathrm{1}}{C_i}}{\displaystyle \frac{\partial C_i}{\partial z}}+({\displaystyle \frac{\mathrm{1}}{H_i}}-{\displaystyle \frac{\mathrm{1}}{H}})+{\mathit{\alpha }}_{T_i}{\displaystyle \frac{\mathrm{1}}{T}}{\displaystyle \frac{\partial T}{\partial z}}].\end{array}$$

Here H is the atmospheric scale height, and ${\mathit{\alpha }}_{T_i}$ is assumed to be zero since the thermal diffusion effect would be of no importance in the stratosphere (Ishidoya et al., 2013). The derived flux was added to the standard transport formulation for each species.

The eddy vertical diffusion in the stratosphere is often neglected in CTMs. However, it should be considered along with molecular diffusion for modeling GS. In the SOCRATES model two options are available for the parameterization of the gravity wave forcing and vertical diffusion coefficient (Huang et al., 1998). The default standard option uses the Lindzen–Holton gravity wave breaking scheme (Lindzen, 1981; Holton, 1982). A second option is to use the parameterization scheme developed primarily by Fritts and Lu (1993), which uses a general gravity wave spectral formulation provided by the observed gravity wave spectrum to deduce how the gravity wave energy flux responds to variation in the background environment. The first scheme was adopted in the current NIES TM simulation. In general, the eddy diffusion mixes concentrations in the volume, reduces vertical stratification, and thereby weakens the molecular diffusion effect, as discussed by Kockarts (2002).

2.3 Simulation setup

The study of GS requires considering two isotopes of atmospheric tracers with a small difference in molecular mass. Following the setup for the SOCRATES baseline-atmosphere run (Ishidoya et al., 2013), we calculated the distributions of ¹³C¹⁶O₂ (M_i=45) and ¹²C¹⁶O₂ (M_i=44). Firstly, a 20-year spin-up calculation with CO₂ using biospheric and oceanic fluxes only and then a 29-year (1988–2016) simulation with total CO₂ fluxes (biospheric, oceanic, and fossil fuel) were performed. The fluxes used were obtained with the GELCA-EOF (Global Eulerian–Lagrangian Coupled Atmospheric model with Empirical Orthogonal Function) inverse modeling scheme (Zhuravlev et al., 2013). This set of fluxes reproduces realistic time and spatial distributions of CO₂ mixing ratio with strong seasonal variations in the Northern Hemisphere (NH) and weak variations in the Southern Hemisphere (SH). The period 1988–2014 is covered by the original JCDAS, which is extended by JRA-55 (an updated version of the Japanese reanalysis) remapped to the same horizontal and vertical grid, as JCDAS production discontinued after 2014. Use of JRA-55 for the whole simulated period is preferable; however, model redevelopment is required to take full advantage of the improved vertical and horizontal resolutions.

The 〈δ〉 value, a measure of the GS, is defined as the isotopic ratio of the CO₂:

$$\begin{array}{}\text{(6)}& \langle \mathit{\delta }\rangle =\mathit{\delta }\left({}^{\mathrm{13}}{\mathrm{C}}^{\mathrm{16}}{\mathrm{O}}_{\mathrm{2}}\right)={\displaystyle \frac{[n\left({}^{\mathrm{13}}{\mathrm{C}}^{\mathrm{16}}{\mathrm{O}}_{\mathrm{2}}\right)/n\left({}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{16}}{\mathrm{O}}_{\mathrm{2}}\right){]}_{\mathrm{strat}}}{[n\left({}^{\mathrm{13}}{\mathrm{C}}^{\mathrm{16}}{\mathrm{O}}_{\mathrm{2}}\right)/n\left({}^{\mathrm{12}}{\mathrm{C}}^{\mathrm{16}}{\mathrm{O}}_{\mathrm{2}}\right)]{}_{\mathrm{trop}}}}-\mathrm{1},\end{array}$$

where subscripts strat and trop denote the stratospheric and tropospheric values, respectively. As the tropospheric value trop we selected the model tracer concentration from the third level, which corresponds to the lower boundary of the free troposphere. The tropospheric 〈δ〉 variations are very small and negligible compared with those in the stratosphere.

CO₂ is a useful tracer of atmospheric dynamics and transport due to its long lifetime. It is chemically inert in the free troposphere and has only a small stratospheric source (up to 1 ppm) from methane oxidation (Boucher et al., 2009). Sufficiently accurate estimation of emissions and sinks, together with knowledge of their trend in combination with the good performance of the NIES model in simulating greenhouse gases, makes the selection of CO₂ appropriate for this study.

3 Results and discussion

3.1 Zonal mean distribution

The zonal mean distribution of CO₂ in the upper part of the atmosphere is driven by the large-scale transport processes: fast quasi-isentropic mixing is combined with upwelling in the tropics and downwelling in the extratropical lowermost stratosphere. In the troposphere, vertical mixing is well developed. With height, the dynamic characteristics weaken, and the mass flux due to molecular diffusion increases (Eq. 3). At a certain level near the tropopause, vertical mixing is no longer able to suppress diffusion and the 〈δ〉 value becomes nonzero.

Figure 1 Annual mean altitude–latitude distributions of 〈δ〉 value (per meg) (a) and the AoA (years) (b) calculated using the NIES TM. The unit (per meg) is typically used in isotope analysis and is the same as 10⁻³ ‰.

Along with the 〈δ〉 value we analyze the AoA (Fig. 1). For this, we used the idealized linearly growing “surface” tracer proposed in the age-of-air intercomparison project (Krol et al., 2018). To fit with our analysis period we extended the original simulation period (1988–2014) to 29 years (1988–2016) with a shorter (10 years) spin-up, as less time is required to reach equilibrium for the AoA analysis.

The processes creating the deformation of the zonal mean cross sections of the GS and the AoA are similar here: the tropical upwelling pumps in tropospheric air and stretches the parameter profile upward (Fig. 1). The NIES TM results shows strong subtropical mixing barriers in both hemispheres compared with the SOCRATES model.

Figure 2 Mean altitude–latitude distributions of the difference between the reciprocals of H_i and H ($L_i^{-\mathrm{1}}$ [1 km⁻¹], see text) for (a) JJA, (b) SON, (c) DJF, and (d) MAM. The results are averaged for 2000–2015.

Figure 3 Mean latitude–longitude distributions of 〈δ〉 value at the 10 hPa level in south polar and north polar projections for (a) JJA, (b) SON, (c) DJF, and (d) MAM. The results are averaged for 2000–2015.

Figure 4 Profiles of 〈δ〉 calculated by the model for February, April, June, August, October, December, and yearly mean averaged for 2000–2015 over the five latitudinal bands shown in Table 1

To estimate the contribution of atmospheric conditions to molecular diffusion, we consider the sum of the three terms in the square bracket of Eq. (5). Because the contribution of the first term (concentration gradient) is relatively small, the second term (originated from pressure gradient in Eq. 1) is the major contributor among the three terms. Therefore, the sum of the three terms can be approximated by the difference between the reciprocals of two scale heights (hereafter referred to as $L_i^{-\mathrm{1}}$). It has a dimension reciprocal to the length and is interpreted as a measure of the efficiency of vertical molecular diffusion under gravity. In view of the essentially one-dimensional nature of GS, it is interesting to consider how $L_i^{-\mathrm{1}}$ distributes in the troposphere and stratosphere. Figure 2 shows the latitude–height distribution of $L_i^{-\mathrm{1}}$ averaged in each season for the case of ¹²C¹⁶O₂. Here the positive values indicate that ¹²C¹⁶O₂ molecules descend relative to major constituents. The temperature fields necessary for the calculation are taken from the JCDAS reanalysis. Since $L_i^{-\mathrm{1}}$ is inversely proportional to temperature, it is generally small in the troposphere and has maxima in cold regions such as the tropical tropopause region and the winter time stratosphere.

The enhancement of $L_i^{-\mathrm{1}}$ does not readily result in a remarkable GS, because it is the difference of $L_i^{-\mathrm{1}}$ between ¹³C¹⁶O₂ and ¹²C¹⁶O₂ that creates GS in our case. We could expect that the enhancement of $L_i^{-\mathrm{1}}$ combined with the long stratospheric transit time in the polar stratosphere will be favorable for GS. Figure 3 compares the horizontal distributions of the seasonal mean 〈δ〉 on 10 hPa pressure surface in polar projections. We can see remarkable GS (small values of 〈δ〉) in both polar regions exhibiting surprisingly clear axial symmetry. In the present analysis, the physical processes that drive GS (Eq. 1) have been rearranged in the form of Eq. (5) to separate the contribution to GS in two factors, one the concentration gradient (the first term) and the other the temperature structure. A stronger seasonal variability of GS in the Southern Hemisphere is caused by changes in vertical pressure gradient (Eq. 1) reflected to those in scale height difference between species.

Table 1 Latitude bands used for averaging 〈δ〉 values.

To minimize local temporal and spatial effects, the seasonal variation of vertical profiles was analyzed for five main climate zones: the southern and northern high latitudes, the southern and northern middle latitudes, and the tropics, as shown in Table 1 and averaged over time (2000–2015). It is clear from Fig. 4 that seasonal variation is evident from a level of about 12 km, except in the tropical region, where it starts from 20 km. The amplitude increases with height and reaches a maximum at the top of the model domain. It is quite obvious that seasonal variability is almost imperceptible in the tropics and increases towards the poles. The maximum and minimum values are reached in December–February and June–August for the South Pole and in April–June and October–December for the North Pole, respectively.

Table 2 Observation sites.

Probably the stronger polar vortex in the SH presumably leads to the enhanced GS (smaller values of 〈δ〉) in JJA in the SH (Figs. 3a and 4a light blue line) compared to that in DJF in the NH (Figs. 3c and 4e). In the NH, on the other hand, GS appears weakest during winter (Fig. 3c) in spite of the enhancement of $L_i^{-\mathrm{1}}$ (Fig. 2c).

Weak sensitivity to seasonal changes of the tracer concentration is a significant advantage of GS over AoA, which is the standard indicator of circulation in the stratosphere. On the other hand, this method requires more accurate sampling tools (i.e., balloon-borne cryogenic samplers) that are more difficult to deploy than other more common methods.

Figure 5 Vertical profiles over Sanriku of 〈δ〉 (per meg) (a, c) and AoA (years) (b, d) calculated from the modeled “surface” tracer (in red) and the observed CO₂ (in blue) distributions for 4 June 2007 (a, b) and 28 August 2000 (c, d).

Figure 6 Same as Fig. 5 but over the Kiruna site for 18 March 1997.

3.2 Vertical profiles

The modeled results are compared with measurements over the main climatic zones: the circumpolar regions, the temperate latitudes, and the tropical latitudes (Table 2). The collection of stratospheric air samples using a balloon-borne cryogenic sampler was initiated in 1985 at the Sanriku Balloon Center of the Institute of Space and Astronautical Science (Nakazawa et al., 1995, 2002; Aoki et al., 2003; Goto et al., 2017; Sugawara et al., 2018; Hasebe et al., 2018). The program has continued up to the present. In addition to observations over Japan, stratospheric air samples were also collected over the Scandinavian Peninsula, Antarctica, and Indonesia.

From those air samples, δ(¹⁵N) of N₂, δ(¹⁸O) of O₂, δ(O₂∕N₂), δ(Ar∕N₂), and δ(⁴⁰Ar) were derived to detect GS in the stratosphere (Ishidoya et al., 2006, 2008, 2013). The effect of GS on the isotopic ratio depends on Δm rather than on the atmospheric component, as follows from Eq. 6. Thus δ values from observations of various tracers and the current simulations can be compared.

Figure 7 Same as Fig. 5 but over the Syowa station for 3 January 1998 (a) and 5 January 2004 (b).

Figure 8 Same as Fig. 5 but over the Biak site for 22–28 February 2015.

Most of the observations were collected in the northern part of Japan over Sanriku (eight profiles) and Taiki (one) in the warm season. Five profiles were observed at the beginning of summer (late May to early June) and four profiles at the end of summer (late August to early September). Typical spring and fall profiles are shown in Fig. 5. For this comparison, the modeled data are daily output at the nearest grid cell.

In the NH, the tropospheric CO₂ is dominated by a strong seasonal cycle due to biospheric activity, which removes CO₂ by photosynthesis during the growing phase to reach a minimum in August–September and releases it during boreal autumn and winter with a maximum in April–May. Due to steady growth and seasonal variation, CO₂ concentrations in the atmosphere contain both monotonically increasing and periodic signals. Spring profiles are smoother, while in autumn they vary with height. The observed AoA shows a strong inversion in the lower altitudes due to seasonal uptake of CO₂, as confirmed by CONTRAIL measurements (Machida et al., 2008; Sawa et al., 2008) and the Lagrangian transport model TRACZILLA (Diallo et al., 2017). The AoA calculated from the modeled surface tracer does not reproduce such effect.

For the high-latitude sites Kiruna and Syowa (Figs. 6–7) the observed profiles are mainly smooth and have smaller vertical fluctuation, apart from the uppermost level over the Syowa station for 5 January 2004.

Due to the limited availability of balloon launching facilities, only one air sample has been conducted in the equatorial mid-stratosphere over Biak (Hasebe et al., 2018). The observed distribution can be explained by the mixing of large-scale NH and SH background values and long-range tracer transport with convective lifting (Inai et al., 2018). For this site, the 〈δ〉 value and the AoA variations with height are very small (Fig. 8), as vertical upwelling pumps young and well-mixed air upward.

Although this is not a model validation paper, it is necessary to evaluate the modeled results by comparison with observations, as the new parameterization for GS simulation was incorporated in the NIES TM. A detailed statistical analysis is summarized for the five observational sites in Table 3. This includes the standard deviation of the misfit between modeled and observed AoA and 〈δ〉 (σ(ΔAoA) and σ(Δ〈δ〉), respectively), the Pearson correlation coefficient between modeled and observed AoA and 〈δ〉 (r(AoA) and r(〈δ〉), respectively), and the Pearson correlation coefficient between AoA and 〈δ〉 from observation (r(AoA, 〈δ〉)_o) and from the model (r(AoA, 〈δ〉)_m). To calculate standard deviation and correlation coefficients only coincident points were selected.

The high values of the correlation coefficients between simulated and observed parameters prove the model efficiency with the implemented parameterization. The lower correlation for GS than for r(AoA) stresses the complexity of the high-precision cryogenic sampling required for GS. For example, most observed profiles show a tendency to have zones of very weak increase or even inversion of the parameters starting from a level of 20–25 km. Despite these limitations, the ability to study the physics underlying GS is a fundamental advantage of the 3-D simulation compared with the 2-D simulation as performed by SOCRATES.

The standard deviations σ(AoA) and σ(〈δ〉) quantify model–observation misfits. We stress a tendency of increase towards the high latitudes, although it seems that a larger gap is obtained for Biak. However, if we normalize the standard deviations by the value of the absolute maximum value of the characteristic, the error decreases towards high latitudes.

Table 3 Standard deviation of misfit between modeled and observed AoA and 〈δ〉 value (σ(ΔAoA) and σ(Δ〈δ〉), respectively); the Pearson correlation coefficient between modeled and observed AoA and 〈δ〉 value (r(AoA) and r(〈δ〉), respectively); the Pearson correlation coefficient between AoA and 〈δ〉 from observation (r(AoA, 〈δ〉)_o) and from the model (r(AoA, 〈δ〉)_m).

GS and AoA are useful indicators of atmospheric transport processes. Two other correlations from Table 3 (r(AoA, 〈δ〉)_o and r(AoA, 〈δ〉)_m) quantify their relationship. Since one value increases with height and the other decreases, the correlation is negative. Modeled results show stronger anti-correlation than that observed, probably due to more straightforward connections in transport simulation; the parameterization used may not take into account additional factors affecting GS. We also do not exclude the influence of observational errors.

3.3 Relationship between age of air and 〈δ〉

Further study of the relationship between the AoA and the 〈δ〉 value is useful for understanding the atmospheric processes, as both would be affected to some extent by changes in the stratospheric circulation. For comparison the modeled AoA and 〈δ〉 values were averaged over the same five broad-latitude bands as in Sect. 3.1 (see Table 1). Along with the modeled values (solid lines) the observed data are also depicted (symbols) in Fig. 9. Observation sites are quite evenly distributed across the selected latitude bands: Syowa in the southern high latitudes, Biak in the tropics, Sanriku and Taiki in the northern middle latitudes, and Kiruna in the northern high latitudes.

Figure 9 shows a near-one-to-one relationship between AoA and GS regardless of latitude and height. Although the model–observation discrepancy is significant for all layers, the distributions of observations show a similar pattern for the lower part of the stratosphere despite large scatter. The positive 〈δ〉 values obtained from observations are not reproducible according to the theory.

The 〈δ〉 value decreases rapidly for age values older than 4 years, as the molecular diffusion coefficient increases with increasing height due to its pressure dependence (Eq. 4), which causes the enhancement of gravitational separation with increasing height. This mechanism does not affect AoA significantly in the stratosphere (Ishidoya et al., 2008, 2013; Sugawara et al., 2018). This emphasizes a nonlinearity in the GS–AoA relationship in the stratosphere.

Another difference stressed the independence of the driving mechanisms of AoA, and GS is the influence of the seasonal variation. For younger air (0–2 years) the seasonal variability of AoA is maximal and falls with altitude, while the 〈δ〉 variability increases continuously upwards, as described in Sect. 3.1.

Sugawara et al. (2018) assumed that vertical upwelling in he tropical tropopause layer (TTL) acts to weaken and downwelling in high latitudes acts to strengthen the effect on the GS. To minimize the discrepancy between the model-calculated (with SOCRATES) and observed 〈δ〉 values over the equatorial region, a correction of the mean meridional circulation and the horizontal eddy diffusivity was performed. However, that correction cannot fully improve the simulation of the stratospheric circulation and hence the reproducibility of AoA and GS by the SOCRATES model. In this study the tracer simulation is based on reanalysis, which describes the atmospheric circulation in sufficient detail. As part of the global balanced air transport the upwelling and downwelling are limited by various constraints including mass flux conservation. The resulting tracer field distribution simultaneously describes AoA and GS, so both parameters can be used to describe the structure of the global atmospheric circulation. However, the limited set of observations and the limitations of the model do not yet allow us to investigate this mechanism and determine the structure of the AoA–GS relationship in more detail.

Chabrillat et al. (2018) presented a consistent intercomparison of AoA according to five modern reanalyses (ERA-Interim, JRA-55, MERRA, MERRA-2, and CFSR) and found significant diversity in the distributions which were obtained with the BASCOE (Belgian Assimilation System for Chemical Observations) transport model, depending on the input reanalysis. They have also found large disagreement between the five reanalyses with respect to the long-term trends of AoA. Thus, an ambitious multi-reanalyses approach is needed to distinguish what is robust in the current GS results from what is not.

Figure 9 Plot of 〈δ〉 value against AoA. The modeled results are averaged for the five latitudinal bands shown in Table 1 and over time for 2000–2015. Symbols represent all available observations for Kiruna, Syowa, Sanriku, Taiki, and Biak sites. Error bar values correspond to 1σ.

4 Conclusions

A three-dimensional simulation of GS in the UTLS zone is performed for the first time using the NIES TM with a molecular diffusion parameterization. We consider the 〈δ〉 values derived from the distribution of two isotopes ¹²C¹⁶O₂ and ¹³C¹⁶O₂. The modeled 〈δ〉 values are compared to observations and the zonal mean distributions from the two-dimensional SOCRATES model.

In comparison with the SOCRATES simulation, the NIES model has a number of significant advantages: a three-dimensional tracer transport simulation driven by global JCDAS reanalysis and a vertical coordinate with isentropic levels. The model is optimized to run the greenhouse gas simulation, as confirmed through various validation and multi-model inter-comparisons. The use of optimized CO₂ fluxes provided realistic tracer distribution and seasonality. Along with these strengths, some weaknesses are also revealed: coarse vertical resolution and the shallow top of the model domain.

The model-to-observation comparison shows that the model with this molecular diffusion parameterization is able to reproduce the mean value and the number of small-scale fluctuations recorded by high-precision cryogenic balloon-borne observations in the lower stratosphere. This reconstruction suggests that the tracer distribution can be explained by the properties of transport, as resolved by meteorological reanalysis and the representation of sub-grid-scale effects as diffusion. Overall, the implemented molecular diffusion parameterization in the NIES TM shows reasonable performance.

We found a strong relationship between the modeled GS and AoA, which is the main indicator of circulation in the stratosphere. However, in contrast to AoA, the GS has a lower sensitivity to seasonal variability, which is a significant issue in studies of atmospheric circulation. Thus, the modeled GS characteristics can provide useful insights and complement AoA information to give a more comprehensive evaluation of structure changes in the UTLS. However, due to the simplified approach and parameterizations, the presented simulation of the GS using the NIES model could not fully achieve the potential of 3-D modeling. Modern reanalysis dataset and recently developed transport models that effectively simulated the upper atmosphere can be employed to address these issues. Since this work is the first in 3-D modeling of GS, we believe this insight is useful for the scientific community working in the field of the UTLS studies.