Stabilization of Isolated Vortices in a Rotating Stratified Fluid

The key element of Geophysical Fluid Dynamics—reorganization of potential vorticity (PV) by nonlinear processes—is studied numerically for isolated vortices in a uniform environment. Many theoretical studies and laboratory experiments suggest that axisymmetric vortices with a Gaussian shape are not able to remain circular owing to the growth of small perturbations in the typical parameter range of abundant long-lived vortices. An example of vortex destabilization and the eventual formation of more intense self-propagating structures is presented using a 3D rotating stratified Boussinesq numerical model. The peak vorticity growth found during the stages of strong elongation and fragmentation is related to the transfer of available potential energy into kinetic energy of vortices. In order to develop a theoretical model of a stable circular vortex with a small Burger number compatible with observations, we suggest a simple stabilizing procedure involving the modification of peripheral PV gradients. The results have important implications for better understanding of real-ocean eddies.


Introduction
Energetic jets and vortices dominate the kinetic energy of flows in planetary atmospheres and the oceans where they are known to play a crucial role in the fluid transport.Understanding their linear and nonlinear stability in a rotating stratified fluid was a subject of many studies in Geophysical Fluid Dynamics.The Gaussian model of pressure anomalies in the radial and vertical directions has been widely used for circular vortices in various theoretical studies as it fits oceanic and laboratory vortices reasonably well.However, such a model has been found to be unstable to small perturbations in the typical parameter range of observed long-lived vortices (see the discussion and references in [1]).The purpose of this work is to develop a smooth stable model of a three-dimensional (3D) axisymmetric vortex in a continuously stratified rotating fluid that can reflect longevity of real-ocean vortices.
The longevity of isolated vortices is closely related to their stability to small perturbations controlled mostly by the distribution of the potential vorticity (PV) along isopycnal surfaces.An isolated circular vortex with a small Rossby number is known to be stable if the radial gradient of PV does not change sign [2].Typically ocean eddies are associated with pronounced PV anomalies (PVA) at the level of maximum azimuthal velocity.The vertical extent of the PVA cores is characterized by a maximum vertical displacement of isopycnal surfaces [3].Baroclinic circular vortices with a constant PVA core, localized both radially and vertically, embedded in uniform PV exterior regions, represent the unconditionally stable states.Such structures have often been used by modelers to study the vortex interactions [4,5].Generally, the presence of opposite PV gradients outside the vortex core does not necessarily imply the growth of perturbations: actual vortex instability depends on the details of the peripheral vortex structure and the ambient stratification.When the Burger number is small (i.e., the eddy is wide compared to the baroclinic radius of deformation) and has no flow at depth (compensated state), most of the theoretical studies with low vertical resolution (e.g., two-layer models) predict the growth of perturbations to an initially circular eddy owing to classical baroclinic instability unless the ocean depth is assumed to be unrealistically large (see [6] and references therein).The addition of a deep co-rotating circulation results in the stabilizing effect in circular eddies for a two-layer ocean model, which can be related to a reduction of opposite PV gradient in the lower layer [7][8][9].Alternatively, the coupling between the upper and lower layers weakens if a third middle layer with uniform PV is added that also reduces the PV gradient in the deep layer at rest and enhances the stability of compensated vortices [1].Here, we explore the more general idea of 3D vortex stabilization for a continuously stratified fluid by correcting peripheral PV gradients.
The nonlinear evolution of unstable circular vortices may result in the formation of various vortex structures including dipoles, tripoles, and multipoles as shown by numerical simulations and laboratory experiments.The possibility of monopolar eddies to transform into multipoles is well known both for two-dimensional (2D) and 3D stratified vortices.The book [5] and recent studies [10][11][12][13] describe various vortical structures and their transformations, and provide a comprehensive account of previous studies in this area.
Substantial progress in high-resolution observations and simulations has revealed an abundance of submesoscale and fine-scale structures in the oceans [14].In particular, a clear evidence of the vertical stacking of sharp variations in density gradients (layering) at the periphery of deep ocean vortices, far from the upper ocean and the bottom, is provided by seismic imaging of the water column (see [15] and references therein).Whereas effects of double-diffusive convection were used for explaining the layering for years [16], recently it has been suggested that isopycnal stirring by the differential rotation in the vortex is able to generate the thin layers seen in seismic imaging [17].On the other hand, stirring can be related to the baroclinic instability of localized eddies with substantial domain of opposite PV gradients below and/or above the vortex core [18,19].
In the present study, we demonstrate the nonlinear 3D evolution of an unstable circular Gaussian anticyclonic vortex in a continuously stratified rotating fluid resulting in vortex elongation and breaking into more intense self-propagating structures.In order to control vortex stability by a reduction of opposite PV gradients, we suggest an iterative method for balanced initialization of a stratified vortex based on prescribing PV along isopycnal surfaces.First, we look at the vortex fragmentation with PV distribution in Gaussian vortices used in previous studies.Second, we modify the initial vortex structure to eliminate opposite isopycnal PV gradients and confirm that such a vortex remains stable for decades in numerical model.Third, we show how the stable eddy can become weakly unstable depending on details of the model setup.
The manuscript is organized as follows.The 3D numerical model and results of simulations for a Gaussian vortex is described in Section 2. The structure of a circular vortex in isopycnal coordinates and an iteration scheme for prescribing the initial vortex for a given PV are presented Section 3. Section 4 presents the conclusions.

Instability of a Gaussian Vortex
To illustrate evolutionary features of vortices related to the details of its radial PV gradient, we present the numerical simulation of an intra-thermocline anticyclonic vortex (meddy) with a Gaussian shape of the angular rotation rate: Fluids 2016, 1, 26 Here, Ro is the Rossby number, f is the Coriolis parameter, r 0 is the vortex radius, h is the vertical scale, r is the horizontal distance from the vortex center, z is the vertical coordinate, and z 0 is the level of maximum azimuthal velocity v = rΩ.The pressure anomaly ρ 0 P (ρ 0 is the reference density) is prescribed from the cyclo-geostrophic balance: The corresponding density field ρ is obtained from the hydrostatic balance gρ where N 0 is the background buoyancy frequency that gives the minimum value of buoyancy frequency N c at the vortex center and the Burger number Bu The Ertel PV in such circular vortex can be written as where the radial gradient of Ω is calculated along isopycnals described by the coordinate ξ (see the next section).The PV anomaly (PVA = PV − f ) is also used further to characterize the vortex structure.
In the following numerical simulation, f = 10 −4 s −1 , r 0 = 30 km, h = 400 m and the meddy is embedded in a uniformly stratified ocean with buoyancy frequency of N 0 = 2.8 × 10 −3 s −1 , so that Bu = 0.14 while Ro = 0.1.For this value of Bu, the growth rate of the fastest growing azimuthal mode (m = 2) in the quasigeostrophic approximation (Ro 1) is about 0.0007 f, according to [18].The numerical model used in this study is the Massachusetts Institute of Technology (MIT) general circulation model [20,21].The computational domain of 480 km × 480 km × 2 km is resolved by the uniform mesh with 256 × 256 × 200 grid points in x, y, and z, respectively.Periodic boundary conditions were assumed in x and y, and the free-slip boundary condition was applied at the bottom.The model was integrated in time with the time-step of ∆t = 100 s.
Radial profiles of initial PV (Equation ( 5)) interpolated into isopycnal surfaces display domains of opposite PV gradients above and below the major vortex core (Figure 1a).The axisymmetric state is an exact steady solution of the model, and the development of small perturbations results from rounding errors in the numerical code.After 300 days of integration, the instability of the Gaussian meddy results in the elongation of the main vortex core of negative vorticity and the splitting of the outer annulus of positive vorticity into two satellites with outward spirals (Figure 2a).Two months later, the elongated core broke into two pieces (Figure 2b), which propagated in opposite directions (Figure 2c,d).This scenario of severe elongation and fragmentation of the vortex core is typical for strongly unstable vortices with dominant azimuthal mode m = 2 in the two-layer models, e.g., [5,8,11,13,22].For continuously stratified fluid, this scenario was described in quasigeostrophic contour-dynamical model as heton instability for vertically offset negative and positive equal vortex cores with uniform PV [23].Less severe instability of a Gaussian vortex may result in shedding smaller vortices with a substantial reduction of the interior circulation of a surviving vortex in 3D simulations [24].
strongly unstable vortices with dominant azimuthal mode m = 2 in the two-layer models, e.g., [5,8,11,13,22].For continuously stratified fluid, this scenario was described in quasigeostrophic contour-dynamical model as heton instability for vertically offset negative and positive equal vortex cores with uniform PV [23].Less severe instability of a Gaussian vortex may result in shedding smaller vortices with a substantial reduction of the interior circulation of a surviving vortex in 3D simulations [24].A remarkable feature of the development of vortex instability is a dramatic growth of the peak vorticity and kinetic energy during the elongation stage and vortex fragmentation (Figure 3).To a lesser degree, the growth of peak vorticity in unstable cyclones was noticed 30 years ago in the two-layer model [22] and for elongating anticyclones in the three-layer model [13].Clearly, it is related to the conversion from the potential energy of the initial vortex into kinetic energy that is typical for baroclinic instability (Figure 3b).Development of sharp vorticity gradients during the elongation stage can be seen around the vortex core in Figure 4a,c which are similar to Figures 5a,c in the three-layer model [13].Figures 4b and 4d illustrate that vortices resulting from the core fragmentation are smaller but rotating more rapidly.Such a mechanism of vortex energizing deserves further investigation.
For small Bu, the major contribution to PV is related to the vortex stretching characterized by ⁄ (the first term in (2.5)).The sign of the vortex stretching changes between the vortex core and upper/lower domains of opposite PVA.For an anticyclonic meddy, the isopycnals in the core are stretched, while in domains above and below the core, they are squeezed (Figure 5).Note that during the elongation stage, domains of vortex squeezing remain vertically aligned with the vortex A remarkable feature of the development of vortex instability is a dramatic growth of the peak vorticity and kinetic energy during the elongation stage and vortex fragmentation (Figure 3).To a lesser degree, the growth of peak vorticity in unstable cyclones was noticed 30 years ago in the two-layer model [22] and for elongating anticyclones in the three-layer model [13].Clearly, it is related to the conversion from the potential energy of the initial vortex into kinetic energy that is typical for baroclinic instability (Figure 3b).Development of sharp vorticity gradients during the elongation stage can be seen around the vortex core in Figure 4a,c which are similar to Figure 5a,c in the three-layer model [13].Figure 4b,d illustrate that vortices resulting from the core fragmentation are smaller but rotating more rapidly.Such a mechanism of vortex energizing deserves further investigation.For small Bu, the major contribution to PV is related to the vortex stretching characterized by N 2 /N 2 0 (the first term in (2.5)).The sign of the vortex stretching changes between the vortex core and upper/lower domains of opposite PVA.For an anticyclonic meddy, the isopycnals in the core are stretched, while in domains above and below the core, they are squeezed (Figure 5).Note that during the elongation stage, domains of vortex squeezing remain vertically aligned with the vortex core (Figure 5a).After fragmentation, formation of layering above smaller vortices is evident in Figure 5b, which is similar to the results presented in [19].Apparently, it arises due to isopycnal stirring by the differential rotation in the vortex as suggested in [17].Note that the vortex stretching is vertically symmetric in the quasigeostrophic model, while in the more accurate Boussinesq model, the upper domain display stronger layering similar to the observation in Figure 2 [19].In this numerical model, the layering diffuses at the later stages while the smaller vortices do not display instability despite the presence of regions with opposite PVA (Figure 5c).Further studies are needed to clarify the stability properties of smaller self-propagating vortices with nearly double values of Ro and Bu.

Theoretical Model of a Stable Vortex
Consider an axisymmetric vortex in an f-plane, ocean model with constant total depth H using the cylindrical isopycnal coordinates (r, ξ), ξ is the unperturbed depth of isopycnal surfaces, which is further used instead of the vertical coordinate.Then, Z(r, ξ) can be introduced for the vertical displacement of isopycnal depths from their unperturbed levels.A typical profile of Z(0, ξ) at the vortex center is shown in Figure 6a (see also Figure 4 in [3]).Then, the Equations ( 2) and ( 3) can be transformed to the form: where = + / + /2 is the anomaly of the Montgomery potential.It has an extremum value ( ) at some isopycnal level that remains horizontal, i.e., Z( , ) = 0.When is near the sea surface, such vortex is surface intensified (which is typical for the ocean rings), while for intrathermocline eddies can be at any depth.The vertical boundary of the vortex core with the main PVA is marked by the isopycnal surface corresponding to maximum Z at the vortex center in the same manner as the radial extent of PVA is marked by the maximum azimuthal velocity ( , ).Here,

Theoretical Model of a Stable Vortex
Consider an axisymmetric vortex in an f -plane, ocean model with constant total depth H using the cylindrical isopycnal coordinates (r, ξ), ξ is the unperturbed depth of isopycnal surfaces, which is further used instead of the vertical coordinate.Then, Z(r, ξ) can be introduced for the vertical displacement of isopycnal depths from their unperturbed levels.A typical profile of Z(0, ξ) at the vortex center is shown in Figure 6a (see also Figure 4 in [3]).Then, the Equations ( 2) and ( 3) can be transformed to the form: where Π = P + ρgz/ρ 0 + N 2 0 z 2 /2 is the anomaly of the Montgomery potential.It has an extremum value Π 0 (ξ 0 ) at some isopycnal level ξ 0 that remains horizontal, i.e., Z(r, ξ 0 ) = 0.When ξ 0 is near the sea surface, such vortex is surface intensified (which is typical for the ocean rings), while for intrathermocline eddies ξ 0 can be at any depth.The vertical boundary of the vortex core with the main PVA is marked by the isopycnal surface corresponding to maximum Z at the vortex center in the same manner as the radial extent of PVA is marked by the maximum azimuthal velocity r c Ω (r c , ξ 0 ).
where = + / + /2 is the anomaly of the Montgomery potential.It has an extremum value ( ) at some isopycnal level that remains horizontal, i.e., Z( , ) = 0.When is near the sea surface, such vortex is surface intensified (which is typical for the ocean rings), while for intrathermocline eddies can be at any depth.The vertical boundary of the vortex core with the main PVA is marked by the isopycnal surface corresponding to maximum Z at the vortex center in the same manner as the radial extent of PVA is marked by the maximum azimuthal velocity ( , ).Here,  The observed structure of ocean rings and intrathermocline eddies can be represented realistically by a vortex with uniform PVA inside the vortex core and zero PVA outside the vortex core [3].In this case, the radial PV gradient does not change sign, which guarantees stability of such a vortex to small perturbations in balanced models.In order to determine the vortex structure for a given PV distribution, we combine Equations ( 7) and ( 9) to get an elliptic equation for Z: which is solved by an iterative method with boundary conditions Z = 0 at the middle level ξ 0 and far from the vortex core r, ξ → ∞ .Corresponding Π and Ω are calculated from Equations ( 6) and (7).
In order to minimize the (potentially destabilizing) influence of the upper and lower boundaries on the eddy, we consider solutions in which isopycnals are horizontal at the sea surface and the bottom.This is accomplished by adding to the solution Equation ( 10) its mirror images with respect to the surface and to the bottom, i.e., symmetric structures centered at −ξ 0 and 2H − ξ 0 .For a Gaussian vortex, the domains with opposite PVA are located below and above the major vortex core (Figure 1a) (see also Figure 2 in [17]), which provide a source of baroclinic instability described in Section 2. This PV pattern was interpolated into isopycnal coordinates and the regions with negative radial PV gradients were effectively eliminated by replacing them with the uniform PV zones (Figure 1b).The resulting solution to Equation (10) has been re-interpolated into the Cartesian coordinate system, and its evolution was simulated using 3D model described in Section 2. While PV-corrected isopycnal displacement does not differ much from that in the Gaussian vortex (Figure 6a), the azimuthal velocity is substantially increased in the PV-corrected case (Figure 6b).Such a solution with uniform PV outside the vortex core corresponds to the two-layer corotating state considered by Benilov [9].Note that without adding mirror images, the isopycnal displacement remains substantial near the surface and the bottom (Figure 6a), leading to radial density gradients near the surface and the bottom, which may also induce Eady-type baroclinic instability.
Figure 7 illustrates that PV-corrected state remains axisymmetric in the numerical simulations for many years.Note that for a smaller computational domain (Figure 8) or without mirror image adjustment (Figure 9), the numerical solutions display weak instability.
Cartesian coordinate system, and its evolution was simulated using 3D model described in Section 2. While PV-corrected isopycnal displacement does not differ much from that in the Gaussian vortex (Figure 6a), the azimuthal velocity is substantially increased in the PV-corrected case (Figure 6b).Such a solution with uniform PV outside the vortex core corresponds to the two-layer corotating state considered by Benilov [9].Note that without adding mirror images, the isopycnal displacement remains substantial near the surface and the bottom (Figure 6a), leading to radial density gradients near the surface and the bottom, which may also induce Eady-type baroclinic instability.
Figure 7 illustrates that PV-corrected state remains axisymmetric in the numerical simulations for many years.Note that for a smaller computational domain (Figure 8) or without mirror image adjustment (Figure 9), the numerical solutions display weak instability.

Figure 1 .
Figure 1.The initial radial PV/f profiles along isopycnals.(a) for a Gaussian vortex (Equation (1)); and (b) for a PV-corrected state described in Section 3.

Figure 2 .
Figure 2. Evolution of a Gaussian meddy.The horizontal sections of relative vorticity at mid-depth (z = z 0 ) normalized by the Coriolis parameter (ς/ f ) are shown at t = 0.9, 1.1, 1.3, and 1.5 years in (a-d), respectively.

FluidsFigure 3 .Figure 4 .
Figure 3.Time records of (a) the maximum relative vorticity normalized by the Coriolis parameter and (b) the net kinetic energy normalized by its initial value.Note the dramatic increase in both maximal vorticity and kinetic energy associated with the eddy breakdown.

Figure 3 .
Figure 3.Time records of (a) the maximum relative vorticity normalized by the Coriolis parameter and (b) the net kinetic energy normalized by its initial value.Note the dramatic increase in both maximal vorticity and kinetic energy associated with the eddy breakdown.

Fluids 2016, 1 Figure 3 .Figure 4 .
Figure 3.Time records of (a) the maximum relative vorticity normalized by the Coriolis parameter and (b) the net kinetic energy normalized by its initial value.Note the dramatic increase in both maximal vorticity and kinetic energy associated with the eddy breakdown.

Figure 4 .
Figure 4. Zonal horizontal sections through the vortex center at t = 1 year (a,c) and t = 2 years (b,d).Relative vorticity at mid-depth (z = z 0 ) normalized by the Coriolis parameter is shown in (a,c).Potential vorticity normalized by the background PV at mid-depth (z = z 0 ) is shown in (b,d).Red curves represent the initial Gaussian vortex.

Fluids 2016, 1 Figure 4 .Figure 5 .
Figure 4. Zonal horizontal sections through the vortex center at t = 1 year (a,c) and t = 2 years (b,d).Relative vorticity at mid-depth ( = ) normalized by the Coriolis parameter is shown in (a,c).Potential vorticity normalized by the background PV at mid-depth ( = ) is shown in (b,d).Red curves represent the initial Gaussian vortex.

Figure 5 .
Figure 5.The vertical sections of vortex stretching N 2 /N 2 0 through the centers of the original vortex and smaller ones after breaking at times corresponding to Figure 2. (a) t = 0.9 years; (b) t = 1.3 years; and (c) t = 2 years.

)Figure 6 .
Figure 6.(a) the isopycnal displacement the vortex axis as a function of depth; and (b) azimuthal velocity at mid-level ( = ) as a function of radius for the Gaussian (red curve), PV-corrected (blue curves) and PV-corrected vortices without mirror image adjustment (green curves).

Figure 6 .
Figure 6.(a) the isopycnal displacement the vortex axis as a function of depth; and (b) azimuthal velocity at mid-level (z = z 0 ) as a function of radius for the Gaussian (red curve), PV-corrected (blue curves) and PV-corrected vortices without mirror image adjustment (green curves).

Figure 7 .
Figure 7. Evolution of the PV-corrected meddy.The horizontal sections of PV at mid-depth ( = ) normalized by the background value (a) initial; (b) t = 15 years; and (c) the vertical section of PV through the vortex center along the dashed line at (b).

Figure 7 .
Figure 7. Evolution of the PV-corrected meddy.The horizontal sections of PV at mid-depth (z = z 0 ) normalized by the background value (a) initial; (b) t = 15 years; and (c) the vertical section of PV through the vortex center along the dashed line at (b).

Fluids 2016, 1 Figure 8 .
Figure 8.The same as in Figure 2 but for the PV-corrected experiment in the small computational domain.Note the weak instability.(a) t = 0.9 years; (b) t = 1.1 years; (c) t = 1.3 years; (d) t = 1.5 years.Figure 8.The same as in Figure 2 but for the PV-corrected experiment in the small computational domain.Note the weak instability.(a) t = 0.9 years; (b) t = 1.1 years; (c) t = 1.3 years; (d) t = 1.5 years.