Estimation of isentropic stirring and mixing and their diagnosis for the stratospheric polar vortex

Abstract. Isentropic stirring and mixing are important processes that determine the distribution of long-lived trace gases in the stratosphere. Stirring stretches tracer contours into filaments and mixing dissipates tracer variance. The combined effects on tracer transport by stirring and mixing are quantified by the effective diffusivity in the modified Lagrangian-mean (MLM) theory that diagnoses tracer transport in an areal coordinate. Here a method is developed to diagnose transport processes based on tracer contours in geographic coordinates. Compared to the MLM theory this method has resolving ability along tracer contours and quantifies stirring and mixing separately. Also, the influence of diabatic motion on the diagnosed stirring is reduced, which is useful for stratospheric analysis where diabatic motion is uncertain. The developed method is validated in a methane simulation experiment. The diagnosed stirring effects are consistent with established knowledge about stratospheric dynamics. Finally, stirring and mixing effects on trace gases in the polar vortex are diagnosed during a northern polar vortex period. According to the diagnosis stirring and mixing always increase the methane concentration in the polar vortex. However, their effects are reversed by vortex movement and deformation in most cases. Only in a few cases, planetary waves can penetrate into the vortex and stirring increases the methane concentration in the vortex. The developed method is readily applicable to diagnose stratospheric transport processes from satellite observed trace gas distributions.



Introduction
Distributions of long-lived trace gases in the stratosphere are determined mostly by the residual circulation and isentropic stirring and mixing (Haynes, 2004). Air following the residual circulation ascends around the tropics, moves polewards and descends in high latitudes. Most long-lived trace gases have chemical sinks in the stratosphere. These chemical sinks decay the tracer and produce gradients in its distribution that reflect the residence time of the air mass in the stratosphere. As a result, the residual circulation always tends to produce a tropics-to-pole gradient in the tracer distribution on each isentropic surface. Along isentropic surfaces, wave disturbances, e.g. planetary waves, stretch and deform the contours of the trace gas. This process is referred to as isentropic diffusion that lead to an irreversible fusion of air masses. Stirring can enhance mixing by producing long, filamentary contours and stronger gradients across the contour. In contrast to the residual circulation, isentropic stirring and mixing tend to decline the tropics-to-pole gradient in the trace gas distributions.
One of the theories that describe these physical processes is the modified Lagrangian-mean (MLM) framework (Nakamura, 1995 and1996). The MLM approach defines an areal/mass coordinate for each value of the mixing ratio as the area/mass enclosed by the corresponding contour. In this way, a functional relationship between the mixing ratio and the areal coordinate can be established. In the two-dimensional case, e.g. along isentropic surfaces, the advection-diffusion equation with respect to spatial coordinates transforms into a one-dimensional diffusion equation with respect to the areal coordinate. In this diffusion equation the combined effect of stirring and mixing is represented by a single parameter referred to as effective diffusivity. It consists of the constant diffusivity of the advection-diffusion equation and a variable scaling referred to as equivalent length. The diffusivity can be estimated by means of the tracer variance equation (Allen and Nakamura, 2001).
The effective diffusivity has been interpreted to diagnose transport barriers and regions with strong mixing in the atmosphere and the ocean (Nakamura and Ma, 1997;Haynes and Shuckburgh, 2000;Allen and Nakamura, 2001;Marshall et al. 2006;Abernathy et al. 2010). One advantage of the MLM theory is that fluxes along the areal coordinate are diagnosed based on the contour of the tracer without using the wind field. This characteristic is valuable in stratospheric applications since accurate wind information are much sparser than trace gas observations in the stratosphere. Further developments of the MLM theory include partitioning of the fluxes into opposing directions (Nakamura, 2004), incorporating diabatic diffusion into the effective diffusivity (Leibensperger and Plumb, 2014) and extending the concept of effective diffusivity to the troposphere (Chen and Plumb, 2014).
A shortcoming of the MLM theory is that the diagnosis based on the areal/mass coordinate does not have resolvability along the tracer contours. This aspect is inconvenient when the theory is applied to trace gases in the stratosphere whose mixing ratios usually present extremes in the tropics. A separate application of the MLM theory to each hemisphere leads to confusion in the tropics since the rising branch of the residual circulation moves North and South during a year. Moreover, trace gases like water vapor and ozone have additional minima in the polar vortex, especially the southern one. In this case, diagnosis defined on the contour of trace gases are averages over well separated regions and not easy to interpret. Attempts to adapt the Lagrangian-mean formalism to a regional mixing diagnostic include Nakamura (2001)  'mixing efficiency' proposed in Nakamura (2001) is the Eulerian-mean counterpart to the equivalent length in the MLM theory. The mixing efficiency relies on a spatial averaging operator over latitudes and longitudes and only represents stirring processes with scales smaller than those of the averaging operator. In contrast, the equivalent length includes all scales above turbulence. The diagnostic by D'Ovidio (2009) combines the tracer-based effective diffusivity and the particlebased Lyapunov exponent calculation and requires more meteorological data than the MLM method.
In this study, we develop a method that diagnoses transport processes with respect to the geographic coordinate based on the tracer contour. The diagnosis has resolvability along the tracer contour and differentiates stirring and mixing. Stirring and mixing have different roles in distributing trace gases in the stratosphere and are controlled by horizontal wind and turbulent processes, respectively. Also, the developed method reduces the influence of diabatic motion on the diagnosed stirring, which is advantageous for stratospheric analysis where observations of diabatic heating are rare.
The remainder of the article is organized as follows: the method is developed and described in Sect.
2. In Sect. 3 we describe the numerical setup that we use as a test case for our method. The analysis of this numerical scenario is described in Sect. 4. A discussion of the method follows in Sect. 5 before we close with a conclusion.

Method
In the stratosphere trace gas contours along isentropic surfaces are modified continuously by horizontal motion along the isentropic surfaces, diabatic motion across the isentropic surfaces, diffusion and chemical sources and sinks. In this section we derive an expression that describes the stretching and deformation of the contour by horizontal motion along isentropic surfaces. We refer to this quantity as (isentropic) stirring. The quantity is based on the evaluation of the contour and does not require isentropic wind information. The derivation of the diagnostic formulas for isentropic stirring is based on the modified Lagrangian-mean framework (e.g. Nakamura, 1995), which we briefly review in the following paragraph.

Modified Lagrangian-mean framework
Consider a long-lived atmospheric trace gas distribution with mixing ratio q=q (x , y ,θ ,t ) , where x and y are geographic coordinates, θ is the vertical potential temperature coordinate and t is time. The temporal evolution of q is described by the advection-diffusion equation where the first term describes advection by isentropic wind along the isentropic surface, the second term describes advection by diabatic (vertical) motion θ and q summarizes all nonconservative processes, particularly turbulent and molecular diffusive mixing, sources and sinks.
Following the concept of the modified Lagrangian-mean (Nakamura, 1995) the density-weighted areal integral of a function f over the the area enclosed by the contour line q on the isentropic surface θ is defined by where σ=−g −1 ∂ p/∂ θ is the potential temperature coordinate density and d(x,y) is an infinitesimal area on the xy-plane. Note that the mass enclosed by the contour line and the isentropic surfaces θ and θ+δ θ is just mδ θ , where m is Generally, the relationship between m and q on an isentropic surface θ at any time t a is one-toone mapping. This property allows to establish the inverse relationship, i.e. q=q (m, θ , t) . As a consequence, functions of q can also be interpreted as functions of m, e.g.
For our derivation the function m(x , y , θ ,t )=m(q ( x , y , θ , t),θ , t) will be of special importance.
It can be interpreted as scaled version of q with the time-varying scaling given by the relation m(q , θ , t) (see Fig. 1 for an illustration).
The modified Lagrangian-mean approach observes the evolution of the atmospheric tracer with respect to this mass coordinate m, rather than in geographic coordinates x and y. The evolution equation of q with respect to m is given by (Nakamura, 1995, Eq. 2.7 (a,b,c)) θ . For a volume bounded by the contour and two isentropic surfaces, the first term on the right hand side represents the diffusive flux across the bounding contour while the other two describe the divergence of diabatic mass flux and vertical transport.
In the following, q is referred to different sets of variables, i.e. (x , y , θ , t) , (q , θ , t) or (m ,θ , t) whenever required. The overline operator, , denotes the global area-weighted average over an isentropic surface.

Derivation of the diagnostic method
Based on the Lagrangian-mean framework, we derive expressions for isentropic stirring and mixing. After a partly technical derivation we present the final forms of both quantities at the end of this section.
The derivation starts by expressing the Eulerian partial time-derivative in Lagrangian-mean coordinates, i.e.
The left hand side of the equation denotes local changes to the trace gas mixing ratio q described by the processes of the Eulerian evolution equation (1). The second term on the right hand side only represents processes that change the distribution of q with respect to the mass coordinate in the Lagrangian-mean evolution equation (2). These nonconservative processes include sources and sinks, diffusive mixing and advection by air flow across isentropic surfaces. Among other processes, isentropic stirring is included in the term ∂m ∂ t | x , y ,θ . To derive an expression for stirring, we rearrange the equation and replace the Eulerian and modified Lagrangian-mean time derivatives with the right hand sides of the evolution equation (1) and (2), respectively. We have The goal is to isolate processes that describe isentropic stirring from those that describe nonconservative and diabatic processes. Recall that we refer to stirring as contour deformation due to motion along isentropic surfaces. The effects by nonconservative processes and diabatic advection should be subtracted from both sides of the equation. After this manipulation there is It is clear that modification to m(x , y , θ ,t ) by nonconservative processes and diabatic advection occurs only when these processes are nonuniform along the contour of q. Eq. 6 is an expression for the contour modification caused by processes other than diffusion, chemical reactions and diabatic advection. Except for vertical advection, diabatic motion can also modify m(x , y , θ ,t ) by changing mean air density over isentropic surfaces. The mean air density variation does not occur explicitly on the right-hand side of Eq. 6. In the following, further analysis is conducted to show that the right-hand side of Eq. 6 consists of horizontal advection along isentropic surfaces and effects by mean air density variations only.
By definition of the operator M and by use of the Leibniz integral rule we have In Eq. 7 ∂σθ ∂θ | x, y ,t contains contributions by horizontal divergent motion caused by vertical convergence of diabatic motion and expansion of air, e.g. due to decreasing air density of uplifting air. Density changes are described by the mass continuity equation. Expressed in the potential temperature coordinate it reads Applying the areal average operator over an isentropic surface, which we denote by an overline, This relation indicates that change in the mean density of an isentopic surface is related to diabatic mass divergence. The horizontal divergent wind ⃗ v d that compensates the local diabatic mass convergence and the expansion of air in an isentropic layer is characterized by In combination with the condition ∇×σ ⃗ v d =0 the horizontal divergent wind is uniquely defined. We use the horizontal divergent wind to reinterpret some terms in Eq. 6. A manipulation of Eq. 6 reads The term defined as D describes advection of m( x , y , θ ,t ) by the horizontal divergent wind where ⃗ n m is the unit vector parallel to the horizontal gradient of m( x , y , θ ,t ) within isentropic surfaces and δ m is a infinitesimal perturbation of m. Sm denotes the area enclosed by the contour line q(m). We use the above relation and subtract the term S m ∂ σ ∂t | θ from both sides of Eq. 9 to receive an expression for isentropic stirring .
The first line represents the diagnostic formula for the stirring effect and the second line explains its physical content. By including the definition of isentropic stirring (Eq. 10) and the evolution equation with respect to m (Eq. 2), Eq. 3 transforms to The right hand side terms in the second line describe the effects of stirring, nonconservative processes, temporal variation of q caused by vertical motion and mean expansion of air, respectively. For a long-lived trace gas with negligible stratospheric sources and sinks, e.g.
methane, molecular and turbulent diffusive mixing are the only nonconservative processes. We neglect the effect of diabatic diffusion and set q=k h ∇ h 2 q with the isentropic diffusivity kh. The effects of isentropic stirring and mixing on the trace gas concentration q are (12) Both diagnostic quantities, isentropic stirring and mixing, are based on the evaluation of the trace gas field q and its scaled version m and do not require isentropic wind information. A procedure to estimate the isentropic diffusivity kh is presented in Section 2.4.

Remarks to isentropic stirring and mixing
The stirring effect, boundary of a stable polar vortex, diffusive effects on m are reduced compared to those on q. The influence of an inaccurately estimated diffusivity is then reduced.

Estimation of isentropic diffusivity
The isentropic and diabatic diffusivities kh and kθ describe the diffusive mixing by molecular and turbulent diffusion along and vertical to the isentropic surfaces. Turbulent diffusion parameterizes the mixing effect from unresolved small scale processes. Global estimates can be inferred from the tracer variance equation (e.g. Leibensperger and Plumb, 2014). Wang et al (2020) extended the method to estimate the average diffusivities for a stratospheric region bounded by two isentropic surfaces θ1 and θ2 by means of the spatially bounded tracer variance equation In this equation M is the airmass of the volume between the two isentropic surfaces θ1 and θ2 and here the overline denotes the mass-weighted average of that volume. kh and kθ are the isentropic and diabatic diffusivities, respectively, and τchem is the chemical lifetime of the trace gas. These

Application
In this section the analysis framework derived in Sect. 2 is applied to estimate stratospheric stirring and mixing from the distribution of a long-lived trace gas in a numerically generated methane scenario. In the following we describe the details of the numerical experiment. the total diabatic heating rates derived from ERA5 temperature tendency properties (Ploeger et al., 2021). The configuration as well as the model initialization follows the model setup described by Pommrich et al. (2014) with 100 km horizontal resolution and 400 m vertical resolution around the tropopause.
The following analysis is restricted to the stratospheric region between 500 K and 1800 K. For the numerical evaluation the modeled daily mean concentration fields are sampled on isentropic surfaces in that region. The estimation of isentropic mixing of the trace gas needs the evaluation of the Laplacian of its mixing ratio within isentropic surfaces. The Laplacian is calculated analytically by expressing the two-dimensional mixing ratio field as a series of spherical harmonics truncated by the triangular T75 truncation (e.g. Daley and Bourassa, 1978). The horizontal resolution of the triangularly truncated approximation is coarser than the model resolution because the effective resolution in simulated mixing ratio fields is typically reduced by diffusion.
We evaluate the stirring from the simulated methane fields using the derived diagnostic formula (Eq. 10a and 12). In that formula, time derivatives are expressed as finite differences of two adjacent daily mean CH4 contours and air density. The other terms require the construction of the

Estimated diffusivities
The isentropic and diabatic diffusivities kh and kθ for the model setup are estimated using the linear regression approach described in Sect. 2.4. Figure 2 present the fitting effects by the forms f(t) = −kθa(t) + c1 and f(t) =−khb(t) + c2. As expected, the distributions of the points (x, y) with f(t) as ycoordinate and a(t) or b(t) as x-coordinate, respectively, show linear relations. This confirms the feasibility of the used method to estimate the diffusivities. The obtained diffusivities for the stratospheric region bounded by 500 K and 1800 K are kh=0.3205×10 5 m 2 /s and kθ=0.0003 K 2 /s.

Reduction of diabatic effects
Contributions by diabatic motion in local temporal variations of q and m are shown in Figure 3. The terms under comparison are θ ∂q ∂θ | x, y ,t and (θ−θ) ∂ q ∂ θ | m ,t , where the local temporal variation of m has been scaled by ∂ q ∂m | θ, t . As the influence of diabatic transport increases with altitude the reduction of such influence is more significant at high levels. For example, diabatic transport of q can be as high as 400 ppb/day in the local temporal variation of q but only around 100 ppb/day in that of m at 1638 K on Jan 27 of 2010. This reduction in diabatic contribution is helpful when deriving horizontal stirring information from satellite-observed trace gases. Due to scarce measurements vertical wind speeds in the stratosphere are highly uncertain and the wind speeds of the modeled residual circulation vary largely between models. The tracer concentration field is determined by 3-dimensional motion of air. The local temporal variation of m is less sensitive to vertical motion than that of q. Therefore, the diagnostic formula for isentropic stirring efficiently extracts horizontal stirring from temporal variations of trace gas concentration created by 3dimensional motion. dominant westerlies in the winter and easterlies in the summer (Andrews et al., 1987).

General characteristics of stirring and mixing
Correspondingly, the rising branch of the residual circulation in the tropics biases towards the summer hemisphere.
CH4 has sources at the surface and undergoes strong oxidation and photolysis reactions in the upper stratosphere and the mesosphere. As a result, upward motion increases and downward motion decreases CH4 mixing ratios in the stratosphere. The rising branch of the residual circulation corresponds to the high-value CH4 mixing ratio region that stretches upward in the tropics (see Fig.   7). The sinking branches indicated by the low-value CH4 mixing ratios stretch downwards from the extratropics towards the polar regions. The sinking branch in hemispheric winter extends to significantly lower altitudes than the sinking branch in hemispheric summer.
The effect of the stirring term on the temporal variation of the methane mixing ratio is calculated using Eqs. 10a and 12. Figure 8 shows the zonal mean distributions of temporal variations of CH4 mixing ratios due to stirring. The dominant effect of stirring is the strengthening of CH4 mixing ratios in sinking branches of the residual circulation. Effects by isentropic stirring in the tropical rising branch are minor. This reflects the fact that the tropical stratosphere is weakly disturbed by extratropical waves due to the subtropical barrier (Chen at al., 1994;Neu et al., 2003). In the extratropics, the most important dynamical processes are the build-up of the polar vortex during winter and its break down during spring. There are lots of wave disturbances during these periods.
The northern polar vortex starts to build up in the middle of SON (Sep., Oct. and Nov.) and breaks down around the end of DJF (Dec., Jan., and Feb.). As a result, stirring effects are large during SON and DJF in the northern extratropics. Similarly, stirring effects are large during JJA (Jun., Jul. and Aug.) and SON in the southern extratropics.
Mixing effects as shown in Fig. 9 are large in regions where strong spatial variations in the trace gas mixing ratios occur, e.g. around transport barriers. In the southern extratropics mixing effects are largest during JJA when the southern polar vortex is stable, i.e. during the period large contrast is present in the trace gas mixing ratios between the inside and the outside of the polar vortex.
Similarly, mixing effects are large during SON in the northern extratropics. The evolution of the vortex area, the mean CH4 mixing ratios in the polar vortex, and the effects of stirring and mixing on the mean CH4 mixing ratios at 781 K are displayed in Fig. 10. The vortex boundary is defined as potential vorticity lines that correspond to the largest potential vorticity gradient with respect to the equivalent latitude. According to the results stirring mostly increases the mean CH4 concentration in the vortex. Mixing increases the mean CH4 concentrations as expected.
However, the combined effects of stirring, mixing and vortex variation mostly decreases the mean CH4 concentration in the vortex. This means that wave disturbances tend to penetrate into the vortex interior but the vortex tends to avoid such penetration through movement and deformation. Also the vortex looses airmasses with high CH4 concentration by adjusting its boundary. In a few cases the penetration into the vortex succeeds and increases the mean CH4 concentration in the vortex.
To explain this point further, Figure 11 shows concentration fields on Jan 29, 2010 (see Fig. 11). On that day a long tail with high-value CH4 mixing ratios penetrates into the vortex. The pattern confirms the diagnosed effect of stirring to increase the methane concentration in the vortex.

Discussion
According to Nakamura (1996) the temporal variation of the trace gas mixing ratio due to stirring and mixing in the two-dimensional case, i.e. within isentropic surfaces, is expressed as where bounded by the actual length of the contour (Haynes and Shuckburgh, 2000). Le is proportional to the length of the contour and the magnitude of the gradient of the trace gas mixing ratio across the contour. The equivalent length is increased by stirring and decreased by diffusion. In this expression, diffusion is included in both the constant and the variable coefficients but stirring is included in the latter only. According to studies by Haynes (2003) andMarshall (2006), dependence of Le on kh varies according to the Pe number. The Pe number describes the relative magnitude of the advective to the diffusive time scale. L e 2 scales like kh -1 when the Pe number is large and approaches a lower limit as the Pe number decreases. As Eq. 14 is expressed with respect to the areal coordinate, the diagnoses of the equivalent length only represents a weighted average along the tracer contour. When applied to methane fields with maxima in the tropics, diagnoses based on the MLM theory is only useful if applied to each hemisphere separately.
Alternatively, a tracer whose concentration decreases/increases monotonically from North to South can be used, e.g. the potential vorticity (Manney and Lawrence, 2016) or an artificial tracer (Allen and Nakamura, 2001).
In the diagnostic quantities developed here, Eqs. 10-12, stirring and mixing are represented by separate terms. The expression for stirring, Eq. 10b, does not depend on the diffusion coefficient kh.
In addition, diagnosed stirring and mixing have resolvability along the contour of the trace gas. This resolvability can be explained as following.
The expression for the stirring effect on m(x , y , θ ,t ) (Eq. 10b) includes two terms, , depends on the local structure of the tracer contour and the magnitude of m that is determined by the nonlocal distribution of the tracer isoline along the whole isentropic surface. After applying the transformation factor ∂ q ∂m | θ, t , as in Eq. 12, the gradient term becomes ∇ q and the nonlocal information is removed. Consequently, the first term of the expression of stirring, is a local term. The second term, D−D , contains the nonuniform part of the divergent motion along the contour and is therefore nonlocal. Its calculation at a particular position is affected by nonlocal processes. In summary, the expression for stirring in Eq. 12 has resolvability along the tracer contour and is local to an extent that the divergent motion induced by vertical mass convergence can be neglected. In Sect. 4 one constant diffusivity is applied to the whole stratospheric region between the isentropic surfaces at 500 K and 1800 K. However, the isentropic diffusivity varies vertically due to vertical variations of wind velocities (Allen and Nakamura, 2001). In their scenario, the estimation of diffusivity is less challenging due to the absence of flow in and out of the region. In addition, diffusivity has inhomgeneities in the horizontal as well as anisotropy (Konopka et al., 2005). Their results reveal that mixing occurs only in regions where flow deforms and shears strongly and diffusivity along the wind is greater than that across the wind. A more precise estimation of diffusivity could be achieved by taking flow deformation into account.

Conclusion
The modified Lagrangian-mean method (