Characteristics of Long-Lived Coherent Vortices in a Simple Model of Quasi-Geostrophic Turbulence †

. Abstract: Macroscale turbulence in the atmosphere is observed to be self-organized into large-scale structures such as zonal jets and robust waves and vortices. A simple model containing the relevant dynamics of turbulence self-organization is quasi-geostrophic turbulence in a stably stratiﬁed atmosphere approximated with a single-layer ﬂuid on a beta-plane. Numerical simulations of this model have shown the dominance of Rossby waves, zonal jets and robust vortices in different regions of the parameter space. In this work, we perform numerical integrations of this model and focus on the regime in which robust large-scale vortices dominate the ﬂow. The goal is to identify the Lagrangian coherent vortices that trap the same air masses in their core throughout their life cycle and to obtain their characteristics. The vortices are identiﬁed using an objective algorithm based on the Lagrangian-averaged vorticity deviation calculated using the advection of Lagrangian particles by the ﬂow. Long-lived vortices with scales comparable to the deformation scale are found with a symmetry between cyclones and anti-cyclones as expected from the simpliﬁed dynamics of the model. The scale as well as the life span of the vortices are also found to increase alongside an increase in the strength of turbulence.


Introduction
Robust, large-scale vortices contribute significantly to the mixing and transport of momentum, heat and constituents in the atmosphere. These vortices also contribute significantly to the observed spatio-temporal variability in the atmospheric circulation, with the stratospheric polar vortex being a prime example [1]. The vortices are supported by atmospheric macro-turbulence through intrinsic processes that are contained even in simplified models of the atmospheric dynamics [2]. One such model utilized to understand the self-organization of turbulence into large-scale structures is the quasi-geostrophic, shallow-water dynamics of a single-layer fluid on a beta-plane with turbulence sustained by random stirring. Numerical simulations of this model have shown that depending on the parameters, the flow is dominated either by zonal jets or by large-scale, isotropic vortices [3]. Previous studies focusing on the vortex regime have provided estimates for the amplitude and the scale of the vortices based on scaling arguments from turbulence cascade theories [4]. However, a thorough study regarding the characteristics of these vortices such as their life span and their propagation properties is lacking, mainly due to difficulties in identifying these structures and following them within the turbulent flow.
In this work, we undertake this task by utilizing a novel technique that objectively identifies Lagrangian coherent vortices, that is, vortices that trap the same air masses in their core throughout their life cycle without any exchanges with their environment [5]. The technique is based on the calculation of the deviation of vorticity along the Lagrangian 2 of 6 trajectories of the air masses. We apply this technique to the turbulent flow produced by the single-layer quasi-geostrophic dynamics in the vortex regime and identify the vortices. We then follow them throughout their life cycle, record their characteristics such as the vorticity at their cores, their scale and their life span, and calculate their statistics.

Quasi-Geostrophic Turbulence in a Shallow-Water, Single-Layer Fluid
Consider a quasi-geostrophic, single-layer fluid on a β-plane. The dynamics can be reduced, in this case, to a single equation for the evolution of potential vorticity [1]: where ψ is the stream function, q = ∇ 2 − λ 2 ψ is the potential vorticity, ∇ 2 = ∂ 2 x +∂ 2 y is the horizontal Laplacian, 1/λ is the Rossby radius of deformation,J (A, B) = ∂ x A∂ y B−∂ y B∂ x A is the Jacobian and β is the gradient of planetary vorticity. Turbulence is sustained by random stirring ξ, representing potential vorticity sources such as convection or processes such as baroclinic instability that are absent from this simplified model. The excitation is assumed to be uncorrelated in time and homogeneous and isotropic in space, injecting energy at a rate ε in a delta ring of radius K f in wavenumber space. To achieve a statistical steady state, there is linear dissipation of the potential vorticity at a rate r.
We integrate (1) in a doubly periodic 2π × 2π channel with a grid of N = 128 points in each direction. The turbulent flow reaches a statistical equilibrium at times 1/r and the integration is carried out until 10/r to ensure stationary statistics. Previous studies showed three different regimes of the turbulent flow that depend on the values of the non-dimensional stratification ∼ λ = λ/K f , the non-dimensional planetary vorticity gradient ∼ β = β/K f r and the non-dimensional strength of the forcing ∼ ε = εK 2 f /r 3 . In the first regime, which is typically found for values of ∼ β ≤ 10, the flow is dominated by large-scale, robust, isotropic vortices [3]. In the other two regimes that are typically found for large values of  Figure 1, where we observe that the flow is populated by an array of vortices. The goal is to identify these vortices using an objective algorithm and to calculate the statistics of their various characteristics such as their vorticity, scale and life span. The technique is based on the calculation of the deviation of vorticity along the Lagrangian trajectories of the air masses. We apply this technique to the turbulent flow produced by the single-layer quasi-geostrophic dynamics in the vortex regime and identify the vortices. We then follow them throughout their life cycle, record their characteristics such as the vorticity at their cores, their scale and their life span, and calculate their statistics.

Quasi-Geostrophic Turbulence in a Shallow-Water, Single-Layer Fluid
Consider a quasi-geostrophic, single-layer fluid on a β-plane. The dynamics can be reduced, in this case, to a single equation for the evolution of potential vorticity [1]: where is the stream function, = (∇ − ) is the potential vorticity, ∇ = is the horizontal Laplacian, 1/ is the Rossby radius of deformation, ( , ) = − is the Jacobian and is the gradient of planetary vorticity. Turbulence is sustained by random stirring , representing potential vorticity sources such as convection or processes such as baroclinic instability that are absent from this simplified model. The excitation is assumed to be uncorrelated in time and homogeneous and isotropic in space, injecting energy at a rate in a delta ring of radius in wavenumber space. To achieve a statistical steady state, there is linear dissipation of the potential vorticity at a rate .
We integrate (1) in a doubly periodic 2 2 channel with a grid of = 128 points in each direction. The turbulent flow reaches a statistical equilibrium at times 1/ and the integration is carried out until 10/ to ensure stationary statistics. Previous studies showed three different regimes of the turbulent flow that depend on the values of the nondimensional stratification = / , the non-dimensional planetary vorticity gradient = / and the non-dimensional strength of the forcing ̃= / . In the first regime, which is typically found for values of 10, the flow is dominated by large-scale, robust, isotropic vortices [3]. In the other two regimes that are typically found for large values of and weak stratification , the flow is dominated by phase-coherent largescale Rossby waves for intermediate values of the excitation strength ̃ [6] and large-scale zonal jets for large values of ̃ [7]. Our focus is on the vortex regime. We therefore choose = 10 and = 4 with = 20 and = 0.1 and integrate the equations for various values of the energy input rate . A typical example of the flow at statistical equilibrium is shown in Figure 1, where we observe that the flow is populated by an array of vortices. The goal is to identify these vortices using an objective algorithm and to calculate the statistics of their various characteristics such as their vorticity, scale and life span.

Vortex Identification Algorithm
To identify the coherent vortices in the flow, we employ the Lagrangian Averaged Vorticity Deviation (LAVD) method developed by Haller et al. [5]. The technique is based on calculating the material contour that nests the same air masses over a specified time interval. This can be achieved by first calculating the Lagrangian paths of air masses that are advected by the flow velocity where x(x 0 , t) is the path of the fluid particle representing an air mass that is initially at position x 0 and u is the flow velocity. We then calculate the LAVD over the specified time interval t 1 -t 0 .
where ζ(x, t) is the vorticity of the fluid particle defined above. LAVD is the vorticity deviation from the domain mean (that is zero for the flow considered) averaged along the particle trajectories and over the specified time interval. The core of the vortices will correspond to the maxima of LAVD, while contours of LAVD surrounding the maxima represent chains of fluid particles around the vortex core that rotate locally at the same rate.
An example of such contours normalized by the square root of spatial mean enstrophy Z m in the flow is shown in Figure 2, along with the fluid particles. The boundary of the vortex is then defined as the outermost convex contour of LAVD. This can be identified by calculating a measure for the dispersion of the particles within a certain closed contour of LAVD. We use a measure that is based on the variance of the positions of the particles within the contour of LAVD and is given by where the bracket denotes the average over all the particles inside the contour. The maximum of this variance within the time interval [t 0 , t 1 ] is then compared to the variance of the same number of particles within a disc of radius R as this is the most compact group of particles. This variance is R 2 /2; therefore, we can define the index measuring their ratio: which is termed the Coherency Index (CI). When the cloud of particles within the LAVD contour is filamented or ejected away from the vortex core, C becomes negative. We therefore set a threshold for the maximum allowed deviation from the dispersion of the particles within the disc. We follow Zhang et al. [8] and set the threshold for C = −3/4, that is, the maximum allowed deviation is set to 75%. The sensitivity tests that were performed showed that the results do not sensitively depend on the exact value for the CI threshold. The coherency index for three contours of LAVD including the identified vortex boundary is also shown in Figure 2. As can be seen, the CI index is almost constant within the vortex core and falls off rapidly outside the vortex boundary.

Vortex Statistics
This objective algorithm is implemented in the following way. Starting from the flow that has reached a statistical equilibrium, we calculate the trajectories of 256 particles for an interval . When the turbulence intensity is large, as in the cases considered in this study, the relevant time scale for the flow evolution is not the dissipation time scale 1/r but the eddy turn-over time scale, which is calculated as = 2 / . As a result, we normalize the time scale by . From these trajectories, the LAVD for the same interval is calculated and its outermost contour, which is the boundary of the vortices, is identified according to the CI threshold. For the array of vortices identified, we calculate their number , the vorticity at their core and their size measured by their diameter . Since some of the vortices are elliptical, an equivalent diameter is calculated for all vortices as = 2 / , where is their surface area. The statistics for the vortices are obtained by slightly perturbing the initial conditions to produce different ensemble members of the flow and repeating the calculations for each member. The number of ensemble members chosen was = 50, as sensitivity studies with more members showed that the results do not change significantly. The calculations are then repeated for larger times . Figure 3 shows the vortex statistics for various time intervals in the case of ̃= 500. Shown are the 25th and 75th percentiles as well as the median of the distribution of the number of the vortices identified to remain coherent until time . Also plotted is the distribution of the vorticity at their cores, which is shown separately for cyclones and anticyclones and the distribution of their diameter. We observe that the number of vortices rapidly decreases as a function of . This means that few of the 17 vortices observed at = 5 remain coherent over much longer time scales. Only about half of them remain coherent at = 20 , only two of them remain coherent at 60 and no vortices can be found for longer time scales. The vorticity at their cores is a few times larger than the average vorticity in the flow and is larger for the few eddies that remain coherent over long times. That is, the long-lived eddies are stronger than their short-lived counterparts. We also observe symmetry between cyclones and anti-cyclones, both in their numbers (not

Vortex Statistics
This objective algorithm is implemented in the following way. Starting from the flow that has reached a statistical equilibrium, we calculate the trajectories of 256 2 particles for an interval τ. When the turbulence intensity is large, as in the cases considered in this study, the relevant time scale for the flow evolution is not the dissipation time scale 1/r but the eddy turn-over time scale, which is calculated as τ e = 2π/ √ Z m . As a result, we normalize the time scale by τ e . From these trajectories, the LAVD for the same interval is calculated and its outermost contour, which is the boundary of the vortices, is identified according to the CI threshold. For the array of vortices identified, we calculate their number N, the vorticity at their core Z and their size measured by their diameter D. Since some of the vortices are elliptical, an equivalent diameter is calculated for all vortices as D = 2 √ A/π, where A is their surface area. The statistics for the vortices are obtained by slightly perturbing the initial conditions to produce different ensemble members of the flow and repeating the calculations for each member. The number of ensemble members chosen was N ens = 50, as sensitivity studies with more members showed that the results do not change significantly. The calculations are then repeated for larger times τ. Figure 3 shows the vortex statistics for various time intervals τ in the case of ∼ ε = 500. Shown are the 25th and 75th percentiles as well as the median of the distribution of the number of the vortices identified to remain coherent until time τ. Also plotted is the distribution of the vorticity at their cores, which is shown separately for cyclones and anti-cyclones and the distribution of their diameter. We observe that the number of vortices rapidly decreases as a function of τ. This means that few of the 17 vortices observed at τ = 5τ e remain coherent over much longer time scales. Only about half of them remain coherent at τ = 20τ e , only two of them remain coherent at 60τ e and no vortices can be found for longer time scales. The vorticity at their cores is a few times larger than the average vorticity in the flow and is larger for the few eddies that remain coherent over long times. That is, the long-lived eddies are stronger than their short-lived counterparts. We also observe symmetry between cyclones and anti-cyclones, both in their numbers (not shown) and in their vorticity. The scale of the vortices roughly coincides with the radius of deformation regardless of the life span of the eddies.
shown) and in their vorticity. The scale of the vortices roughly coincides with the radius of deformation regardless of the life span of the eddies.  Figure 4 shows the vortex statistics for stronger excitation ̃= 5000. While the number of vortices again decreases with , the drop is less rapid than in the case of weaker excitation, and there are a few vortices with longer lifetimes. The size of the eddies remains unchanged and equal to the deformation scale, but the vorticity at their cores is slightly larger in the case of stronger excitation, especially for the long-lived vortices. In summary, we have longer-lived and slightly stronger vortices for a larger turbulence intensity, a result that holds for larger ̃ as well (not shown).   Figure 4 shows the vortex statistics for stronger excitation ∼ ε = 5000. While the number of vortices again decreases with τ, the drop is less rapid than in the case of weaker excitation, and there are a few vortices with longer lifetimes. The size of the eddies remains unchanged and equal to the deformation scale, but the vorticity at their cores is slightly larger in the case of stronger excitation, especially for the long-lived vortices. In summary, we have longer-lived and slightly stronger vortices for a larger turbulence intensity, a result that holds for larger   Figure 4 shows the vortex statistics for stronger excitation ̃= 5000. While the number of vortices again decreases with , the drop is less rapid than in the case of weaker excitation, and there are a few vortices with longer lifetimes. The size of the eddies remains unchanged and equal to the deformation scale, but the vorticity at their cores is slightly larger in the case of stronger excitation, especially for the long-lived vortices. In summary, we have longer-lived and slightly stronger vortices for a larger turbulence intensity, a result that holds for larger ̃ as well (not shown).

Conclusions
In this work, a simple model of quasi-geostrophic turbulence was analyzed with the goal of obtaining the statistics of robust, large-scale, coherent vortices that dominate the flow in a certain parameter regime. The Lagrangian coherent vortices, that is, the vortices keeping the same air masses trapped within their cores, were objectively identified using the Lagrangian Averaged Vorticity Deviation method. Using this technique, a vortex is identified in terms of the vorticity deviation of Lagrangian particles advected by the flow over a specified time interval. The maximum of the vorticity deviation represents the vortex core and the boundary of the vortex is calculated as the outermost convex contour surrounding the maximum. This contour is identified using a threshold related to the maximum allowed dispersion of particles inside the vortex boundary relative to the dispersion of the most compact cloud of particles. Using an algorithm that applies this technique and integrating the equations with slightly different initial conditions to obtain several independent ensemble members, we calculated the average number of vortices, their average size and their average vorticity at their cores for various levels of turbulence intensity. The flow was found to be populated by a number of vortices. Only a few of them have long life spans, and these long-lived vortices appear to be stronger than their counterparts with shorter life spans. Symmetry in the number of cyclones and anti-cyclones was found to hold (at least in this simple model of geostrophic turbulence) with the two species also having the same characteristics. The average size of the vortices was found to be comparable to the Rossby radius of deformation. In addition, the vortices became longer-lived and slightly stronger when the turbulence intensity was larger.
While the simple model used in this work does not target a realistic simulation of atmospheric circulation, the turbulent dynamics maintaining the coherent vortices are universal and inherent even in realistic flows. Therefore, the results of this study can shed light on the inner workings of the stratospheric polar vortex, which exerts great influence on atmospheric circulation [9] but also traps atmospheric constituents, affecting the local chemistry and the radiation budget in the stratosphere [10].