Stable Vortices in a Continuously Stratified Ocean with Thin Active Layer

Abstract

This paper presents a model which yields examples of stable vortices in a continuously stratified rotating fluid, thus providing a possible explanation of the observed longevity of oceanic eddies. The model is based on two assumptions. Firstly, the ocean comprises a thin upper (active) layer and a thick lower (passive) one, with large and small vertical gradients of density, respectively. Secondly, the Rossby number is small, justifying the use of the geostrophic and quasi-geostrophic approximations for the active and passive layers (the two are treated differently because the vortex-induced displacement of the isopycnal surfaces is comparable to the depth of the active layer, but is much smaller than that of the passive one). Using the asymptotic equations derived on the basis of the above assumptions, we prove a stability criterion and thus identify a class of stable vortex profiles. This class is much wider than the one following from the standard requirement that the potential vorticity be monotonic in the whole bulk of the fluid.

1. Introduction

It has been known since the 1977 paper by Lai and Richardson [1] that oceanic eddies can exist for years, but most of the subsequent theoretical work [2,3,4,5,6,7,8,9,10,11] seems to suggest that vortices with "real" parameters should be unstable. The instability is caused by the large radii of most observed eddies (making their Burger numbers $Bu$ small), as large vortices—with some exceptions listed below—tend to be unstable.

The first step toward resolving the paradox was made in 1995 by Dewar and Killworth [12]. They examined a Gaussian vortex in the upper layer of a two-layer ocean and a co-rotating flow in the lower layer: it turned out that a relatively weak “deep flow” could stabilize the vortex. This effect has been explored for other vortex profiles in [13], but none except the Gaussian becomes entirely stable.

The case in which the deep circulation has the same radial profile as the upper-layer vortex was considered in [12,13]. A different approach was used in [14,15]: these papers examined the deep flow that makes the lower-layer potential vorticity (PV) uniform and showed that it stabilizes all vortices, not just the Gaussian one. It was also argued in [14,15] that “deep flow with uniform PV arises naturally below oceanic eddies, as suggested by the fact that most of them are shed by unstable frontal currents. Accordingly, when a vortex moves away from the current and arrives to a new location, the PV field below it cannot change and remains equal to its initial, background value”.

The main shortcoming of the four papers [12,13,14,15] in which stable vortices were found is that they all used the (relatively crude) two-layer approximation of the ocean. In order to model real oceanic eddies, one should instead use a continuous approximation.

This problem was first addressed in [16], where it was shown that a vortex with a constant PV value, surrounded by fluid with a different (but also constant) PV value, is stable. This result fits nicely with the concept of “uniform PV in the lower layer” advocated in [14,15], and, thus, appears to resolve the paradox of the longevity of oceanic eddies. On the other hand, non-uniform-PV vortices with continuous stratification turned out to be unstable [17,18,19].

However, in [16], the stability of uniform-PV vortices in a uniform-PV fluid was proved, but only for the latter being unbounded. The proof does not seem to have an obvious extension to the case of a finite-depth layer, which would be a more realistic model of the ocean.

In fact, there might be a good reason for the lack of a stability criterion for finite layers, as at least some of the vortices considered in [16]—if placed between rigid boundaries—become unstable. As shown in [8], all vortices with $Bu\ll 1$ (including those with uniform PV, in a uniform-PV fluid) are baroclinically unstable, subject to two conditions: (1) if they involve a sizeable vertical shear; and (2) if their thickness is of the order of the ocean’s depth. If at least one of these requirements is violated, the growth rate of the instability becomes small, and one has to take into account the effect of critical levels (which were neglected in [8]). Thus, given that real oceanic eddies are typically strongly sheared, they can be stable only if they are thin—and will be even more so, as the effect of critical levels on the stability of thin vortices is indeed strong [9,10,11].

The present paper investigates the stability of thin vortices in a continuously stratified fluid bounded by two horizontal planes. In Section 2, we formulate the problem mathematically and, in Section 3, we derive an asymptotic model for the limit described in the abstract. In Section 4, the model is used to prove a stability criterion and, thus, identify a wide class of stable thin vortices.

2. The Formulation

We consider a layer of an ideal fluid characterized by the density $\rho$, pressure p, and velocity components $\left(u,v,w\right)$—all depending on the Cartesian coordinates $\left(x,y,z\right)$ and time t. We assume that the fluid is bounded by two horizontal rigid planes located at $z=0$ (bottom) and $z=H$ (surface), and that it is rotating about a vertical axis with an angular velocity f (in the oceanic context, $f=2\Omega \cos \lambda$, where $\Omega$ is the angular velocity of the earth’s rotation and $\lambda$ is the local latitude). We also introduce the fluid’s mean density ${\rho }_0$, and the acceleration due to gravity g.

2.1. The Governing Equations

We employ the Euler equations under the hydrostatic approximation, the f-plane approximation of rotation, and the Boussinesq approximation of the density stratification:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}+\frac{1}{{\rho }_0}\frac{\partial p}{\partial x}=fv$$

(1)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}+\frac{1}{{\rho }_0}\frac{\partial p}{\partial y}=-fu$$

(2)

$$\frac{\partial p}{\partial z}=-g\rho$$

(3)

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$$

(4)

$$\frac{\partial \rho }{\partial t}+u\frac{\partial \rho }{\partial x}+v\frac{\partial \rho }{\partial y}+w\frac{\partial \rho }{\partial z}=0$$

(5)

Equations (1)–(5) are to be solved with the no-through-flow condition at the rigid boundaries:

$$w=0\mathrm{at}z=0,H.$$

(6)

In what follows, in is convenient to use the "semi-Lagrangian" coordinates, such that the vertical coordinate z is replaced by a Lagrangian marker, whereas x and y remain Eulerian. Such coordinates have been used previously in [20,21,22] for shallow-water models similar to ours, and in [23], for a model comprising ours as a particular case. Similar coordinates are often used for general oceanographic modeling (e.g., [24]), where they are referred to as “isopycnic”.

The semi-Lagrangian coordinates are denoted by $\left(t_{new},x_{new},y_{new},\xi \right)$, and we let them be related to the Eulerian coordinates by

$$t=t_{new},x=x_{new},y=y_{new},z=Z(t_{new},x_{new},y_{new},\xi )$$

where the function Z satisfies the following (the subscript ${}_{new}$ omitted):

$$\frac{\partial Z}{\partial t}+u\frac{\partial Z}{\partial x}+v\frac{\partial Z}{\partial y}-w=0$$

(7)

$$Z=0\mathrm{at}\xi =0$$

(8)

$$Z=H\mathrm{at}\xi =H$$

(9)

Rewriting Equations (1)–(6) in terms of the new variables, and taking into account Equation (7), we obtain the following (${}_{new}$ omitted):

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\frac{1}{{\rho }_0}\left(\frac{\partial p}{\partial x}-\frac{\partial Z}{\partial x}\frac{1}{n}\frac{\partial p}{\partial \xi }\right)=fv$$

(10)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\frac{1}{{\rho }_0}\left(\frac{\partial p}{\partial y}-\frac{\partial Z}{\partial y}\frac{1}{n}\frac{\partial p}{\partial \xi }\right)=-fu$$

(11)

$$\frac{\partial p}{\partial \xi }=-\rho gn$$

(12)

$$\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)n-\frac{\partial Z}{\partial x}\frac{\partial u}{\partial \xi }-\frac{\partial Z}{\partial y}\frac{\partial v}{\partial \xi }+\frac{\partial w}{\partial \xi }=0$$

(13)

$$\frac{\partial \rho }{\partial t}+u\frac{\partial \rho }{\partial x}+v\frac{\partial \rho }{\partial y}=0$$

(14)

where

$$n=\frac{\partial Z}{\partial \xi }$$

(15)

Equation (13) can be rewritten in a more convenient form if one adds to it the derivative of Equation (7) with respect to $\xi$, which yields

$$\frac{\partial n}{\partial t}+\frac{\partial \left(un\right)}{\partial x}+\frac{\partial \left(vn\right)}{\partial y}=0$$

(16)

Equation (14) can, in turn, be satisfied by letting

$$\rho =\rho \left(\xi \right).$$

Physically, $\rho \left(\xi \right)$ describes the density stratification of the unperturbed ocean and, thus, can be treated as a known parameter. We also assume that the stratification is statically stable, that is,

$$\frac{\mathrm{d}\rho }{\mathrm{d}\xi }<0.$$

With $\rho \left(\xi \right)$ being the unperturbed stratification and Z describing the depths of isopycnal surfaces in a perturbed ocean, $Z-\xi$ represents the displacement of isopycnal surfaces.

Next, we substitute $\partial p/\partial \xi$ from Equation (12) into Equations (10) and (11), which yields

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+\frac{\partial \psi }{\partial x}=fv$$

(17)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+\frac{\partial \psi }{\partial y}=-fu$$

(18)

where the modified pressure variable can be defined as

$$\psi =\frac{p}{{\rho }_0}+\frac{g}{{\rho }_0}\left[\rho z-{\int }_{\xi }^H{\xi }^{\prime }\frac{\mathrm{d}\rho \left({\xi }^{\prime }\right)}{\mathrm{d}{\xi }^{\prime }}\mathrm{d}{\xi }^{\prime }\right]$$

(19)

Note that $\psi$ involves a degree of arbitrariness, as adding to it a function of $\xi$ and t does not change Equations (17) or (18).

To relate $\psi$ and n, we differentiate Equation (19) with respect to $\xi$ and, recalling Equations (12) and (15), obtain

$$\frac{\partial \psi }{\partial \xi }=-\left(Z-\xi \right)N^2$$

(20)

where

$$N\left(\xi \right)={\left(-\frac{g}{{\rho }_0}\frac{\mathrm{d}\rho }{\mathrm{d}\xi }\right)}^{1/2}$$

is the Väisälä–Brunt frequency. As follows from Equations (20) and (8)–(9),

$$\frac{\partial \psi }{\partial \xi }=0\mathrm{at}\xi =0,H$$

(21)

Finally, we divide Equation (20) by $N^2$ and differentiate with respect to $\xi$ again, which yields

$$n=1-\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }\right)$$

(22)

Equations (16)–(18), (21) and (22) form a governing set for n, u, v, and $\psi$.

We also need an expression for the displacement of fluid particles, $Z-\xi$, in terms of $\psi$. It follows from Equation (20) that

$$Z-\xi =-\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }.$$

In what follows, we examine situations involving a discontinuity (jump) in $N\left(\xi \right)$. To extend our model to such cases, we require that the pressure and the particles’ displacements be continuous across the jump. Mathematically, these requirements amount to

$${\left(\psi \right)}_{\xi ={\xi }_j-0}={\left(\psi \right)}_{\xi ={\xi }_j+0},{\left(\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }\right)}_{\xi ={\xi }_j-0}={\left(\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }\right)}_{\xi ={\xi }_j+0}$$

(23)

where ${\xi }_j$ is the jump’s location.

2.2. Discussion

(1) The mere fact that $\rho$ depends only on $\xi$ means that, at the boundaries ($\xi =0,H$), the density assumes constant values. This restriction appears to rule out modeling near-surface fronts and eddies, as the initial condition for $\rho$ in this case does depend on x, y, and t. This problem can be resolved in two different ways. Firstly, one can let $\rho$ depend on all four variables and, thus, retain Equation (14) in the governing set. Secondly, one can introduce infinitesimally thin layers near the boundaries, within which the density rapidly changes from the value prescribed by the initial condition to a certain fixed value (for the upper/lower boundary, this value should be the minimum/maximum of $\rho$ in the initial condition). Because the "added" layers have zero mass, they do not affect the global dynamics, but still make $\rho$ at the boundaries constant.

The fact that we have the set $\rho =\rho \left(\xi \right)$ implies that we have chosen the latter approach.

(2) The change from the Eulerian to semi-Lagrangian variables is non-singular only if the isopycnal surfaces never overturn (fold)—otherwise, there would be multiple values of Z for the same $\xi$. This restriction is implied everywhere in this work.

3. Asymptotic Equations

3.1. The Assumptions

We let the ocean comprise two layers, of depths $H_1$ and $H_2$, such that

$$\frac{H_1}{H}=\epsilon ,\frac{H_2}{H}=1-\epsilon$$

(24)

where $\epsilon \ll 1$. Thus, the upper (active) layer is assumed to be thin.

Despite the depth difference, we let the two layers involve comparable density variations. Such a situation is commonplace in the real ocean, where a relatively thin upper layer (∼400–$700 \mathrm{m}$) comprises as much density variation as the rest of the ocean. The same applies to near-surface eddies and fronts, as they are stratified differently from the well-mixed water masses beneath them, and rarely penetrate deeper than $500 \mathrm{m}$. In all such cases, one can assume that the upper-layer Väisälä–Brunt frequency N exceeds the lower-layer N by an order of magnitude. Furthermore, if the transitional region between the layers is thin, it can be approximated by a jump in $N\left(\xi \right)$, located at ${\xi }_j=H_2$.

As for the flow’s horizontal scale L and velocity scale U, we assume these to be global, i.e., independent of the layer.

Using the parameters introduced, one can define the Rossby number (for both layers) and two Burger numbers:

$$Ro=\frac{U}{fL},B{u}_1=\frac{H_{1}g\Delta \rho }{f^2{L}^2{\rho }_0},B{u}_2=\frac{H_{2}g\Delta \rho }{f^2{L}^2{\rho }_0}$$

(25)

where the density variation $\Delta \rho$ can be chosen as, say,

$$\Delta \rho =\rho \left(0\right)-\rho \left(H\right)$$

(26)

Evidently, $B{u}_1=\epsilon B{u}_2/\left(1-\epsilon \right)$, and hence, the Burger numbers cannot both be of order-one. We assume that $B{u}_1\ll 1$, which corresponds to the most interesting kind of oceanic eddies—the large eddies—which appear to be theoretically unstable, but sometimes persist for years. No assumption is made regarding $B{u}_2$.

In a multi-parameter asymptotic problem such as ours, it is convenient to inter-relate the parameters, targeting one of the available “characteristic limits”. In the present case, there seems to be only one of these, namely,

$$Ro=\epsilon ,B{u}_1=\epsilon ,B{u}_2=1-\epsilon .$$

Given the equalities of Equation (25), it is implied that

$$L=\frac{1}{f}{\left(\frac{Hg\Delta \rho }{{\rho }_0}\right)}^{1/2},U=H_1{\left(\frac{g\Delta \rho }{H{\rho }_0}\right)}^{1/2}$$

(27)

The use of a single velocity scale for both layers excludes the flows that are fully localized in the upper layer. The flows included must to some extent penetrate the lower layer (but may of course decay toward the ocean’s bottom).

Thus, the separation of the ocean into an active layer and a passive layer is made on the basis of the density stratification, rather than the intensity of the flow.

3.2. The Derivation

According to the assumptions above, the Lagrangian marker $\xi$ in the range $\left(0,H_2\right)$ should be scaled by the passive layer’s depth $H_2$, and in the range $\left(H_2,H\right)$, by $H_1$. Marking the scaled $\xi$ (and other scaled variables to be introduced later) with subscripts ${}_s$, we have

$${\xi }_s=\left\{\begin{array}{cc}{\displaystyle \frac{\xi -H_2}{H_2}}\hfill & \mathrm{if}\xi \in \left(0,H_2\right),\\ \hfill {\displaystyle \frac{\xi -H_2}{H_1}}& \mathrm{if}\xi \in \left(H_2,H\right).\end{array}\right.$$

As the Väisälä–Brunt frequency involves a derivative with respect to $\xi$, it should also be scaled depending on the layer:

$$N_s=\left\{\begin{array}{cc}{\left({\displaystyle \frac{H_2{\rho }_0}{g\Delta \rho }}\right)}^{1/2}N\hfill & \mathrm{if}\xi \in \left(0,H_2\right)\\ \hfill {\left({\displaystyle \frac{H_1{\rho }_0}{g\Delta \rho }}\right)}^{1/2}N& \mathrm{if}\xi \in \left(H_2,H\right)\end{array}\right.$$

where $\Delta \rho$ is defined by (26). The rest of the variables are scaled as follows:

$$t_s=\frac{U}{L}t,\left(x_s,y_s\right)=\frac{\left(x,y\right)}{L},$$

$$n_s=n,\left(u_s,v_s\right)=\frac{\left(u,v\right)}{U},{\psi }_s=\frac{\psi }{fUL},$$

where L and U are defined by Equation (27).

In terms of the scaled variables, Equations (16)–(18) and (21)–(23) with ${\xi }_j=H_2$ become (the subscript ${}_s$ omitted):

$$\frac{\partial n}{\partial t}+\frac{\partial \left(un\right)}{\partial x}+\frac{\partial \left(vn\right)}{\partial y}=0$$

(28)

$$\epsilon \left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)+\frac{\partial \psi }{\partial x}=v$$

(29)

$$\epsilon \left(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)+\frac{\partial \psi }{\partial y}=-u$$

(30)

$$n=\left\{\begin{array}{cc}\hfill 1-{\displaystyle \frac{\epsilon }{1-\epsilon }}{\displaystyle \frac{\partial }{\partial \xi }}\left({\displaystyle \frac{1}{N^2}}{\displaystyle \frac{\partial \psi }{\partial \xi }}\right)& \mathrm{if}\xi \in \left(-1,0\right)\\ 1-{\displaystyle \frac{\partial }{\partial \xi }}\left({\displaystyle \frac{1}{N^2}}{\displaystyle \frac{\partial \psi }{\partial \xi }}\right)\hfill & \mathrm{if}\xi \in \left(0,1\right)\end{array}\right.$$

(31)

$$\frac{\partial \psi }{\partial \xi }=0\mathrm{at}\xi =\pm 1$$

(32)

$${\left(\psi \right)}_{\xi =-0}={\left(\psi \right)}_{\xi =+0},{\left(\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }\right)}_{\xi =-0}={\left(\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }\right)}_{\xi =+0}$$

(33)

where $\epsilon$ is defined by Equation (24), and the ranges $\xi \in \left(-1,0\right)$ and $\xi \in \left(0,1\right)$ correspond to the lower and upper layer, respectively.

Following the usual routine of deriving the (quasi)geostrophic approximation, consider the following combination of the governing equations:

$$\frac{\partial }{\partial y}(29)-\frac{\partial }{\partial x}(30)+(28),$$

which yields the vorticity equation

$$\epsilon \left[\frac{\partial }{\partial y}\left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}\right)-\frac{\partial }{\partial x}\left(\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}\right)\right]+\frac{\partial n}{\partial t}+\frac{\partial \left(un\right)}{\partial x}+\frac{\partial \left(vn\right)}{\partial y}=\frac{\partial v}{\partial y}+\frac{\partial u}{\partial x}$$

(34)

Equation (34) will be used instead of Equation (28).

Next, we let

$$n=\left\{\begin{array}{cc}1+\epsilon \eta \hfill & \mathrm{if}\xi \in \left(-1,0\right)\\ 1+\eta \hfill & \mathrm{if}\xi \in \left(0,1\right)\end{array}\right.$$

Now, keeping only the leading-order terms of Equations (29)–(34), we obtain

$$v=\frac{\partial \psi }{\partial x},u=-\frac{\partial \psi }{\partial y}$$

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial }{\partial t}}\left(\eta -{\nabla }^2\psi \right)+J\left(\psi ,\eta -{\nabla }^2\psi \right)& =0\mathrm{if}\xi \in \left(-1,0\right)\phantom{{\displaystyle \underset{a}{\int }}}\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \eta }{\partial t}}+J\left(\psi ,\eta \right)& =0\mathrm{if}\xi \in \left(0,1\right)\hfill \end{array}$$

(36)

$$\eta =-{\displaystyle \frac{\partial }{\partial \xi }}\left({\displaystyle \frac{1}{N^2}}{\displaystyle \frac{\partial \psi }{\partial \xi }}\right)$$

(37)

$$\frac{\partial \psi }{\partial \xi }=0\mathrm{at}\xi =\pm 1$$

(38)

$${\left(\psi \right)}_{\xi =-0}={\left(\psi \right)}_{\xi =+0},{\left(\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }\right)}_{\xi =-0}={\left(\frac{1}{N^2}\frac{\partial \psi }{\partial \xi }\right)}_{\xi =+0}$$

(39)

where the "horizontal" Laplacian and Jacobian are

$${\nabla }^2\psi =\frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2},J(\psi ,\eta )=\frac{\partial \psi }{\partial x}\frac{\partial \eta }{\partial y}-\frac{\partial \psi }{\partial y}\frac{\partial \eta }{\partial x}$$

In what follows, the asymptotic Equations (35)–(39) will be used for examining the stability of vortices.

3.3. Discussion

Equations (35)–(38) resemble the standard QG set, but—as far as the upper layer is concerned—the resemblance is superficial. Our upper-layer Equation (36) accounts for the order-one displacement of isopycnal surfaces, whereas its QG counterpart is based on the assumption that it is small.

Moreover, even if the usual QG description is relaxed (so that the upper-layer isopycnal displacement can be large), the resulting model would involve a moving interface between the two layers, complicating the problem mathematically.

In the present model, this issue has been resolved by using the semi-Lagrangian coordinates, so that Equations (35) and (36) are solved in fixed domains.

4. Stability of Vortices

4.1. The Linearized Problem

When dealing with axisymmetric vortices, it is convenient to use the polar coordinates $\left(r,\theta \right)$ related to the Cartesian coordinates by

$$x=r\cos \theta ,y=r\sin \theta .$$

Vortices are described by the following solution of Equations (35)–(39):

$$\psi =\mathrm{\Psi }(r,\xi )$$

where $\mathrm{\Psi }$ satisfies

$$\frac{\partial \mathrm{\Psi }}{\partial \xi }=0\mathrm{at}\xi =\pm 1$$

and

$${\left(\mathrm{\Psi }\right)}_{\xi =-0}={\left(\mathrm{\Psi }\right)}_{\xi =+0},{\left(\frac{1}{N^2}\frac{\partial \mathrm{\Psi }}{\partial \xi }\right)}_{\xi =-0}={\left(\frac{1}{N^2}\frac{\partial \mathrm{\Psi }}{\partial \xi }\right)}_{\xi =+0}$$

(40)

It is also implied that

$$\begin{array}{cc}\hfill \frac{\partial \mathrm{\Psi }}{\partial r}=\mathcal{O}\left(r\right)& \mathrm{as}r\to 0,\phantom{{\displaystyle \underset{A}{\int }}}\hfill \\ \hfill \frac{\partial \mathrm{\Psi }}{\partial r}\to 0& \mathrm{as}r\to \infty ,\hfill \end{array}$$

where the latter condition is self-explanatory, and the former guarantees that the fluid’s angular velocity is finite at the center of the vortex:

$$\Omega =\frac{1}{r}\frac{\partial \mathrm{\Psi }}{\partial r}$$

(41)

To examine the stability of the vortex solution, we let

$$\psi =\mathrm{\Psi }+\tilde{\psi }(t,r,\theta ,\xi ),\eta =-\frac{\partial }{\partial \xi }\left[\frac{1}{N^2}\frac{\partial \left(\mathrm{\Psi }+\tilde{\psi }\right)}{\partial \xi }\right]$$

(42)

where $\tilde{\psi }$ describes a small perturbation. Substituting Equation (42) into Equations (35)–(39) and linearizing, we obtain

$$\begin{array}{c}\left({\displaystyle \frac{\partial }{\partial t}}+\Omega \frac{\partial }{\partial \theta }\right)\left[\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \tilde{\psi }}{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \tilde{\psi }}{\partial r}\right)+\frac{1}{r^2}\frac{{\partial }^2\tilde{\psi }}{\partial {\theta }^2}\right]\hfill \\ -\frac{\partial \tilde{\psi }}{\partial \theta }\left\{\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial \left(r^2\Omega \right)}{\partial r}\right]\right\}=0\mathrm{if}\xi \in \left(-1,0\right)\hfill \end{array}$$

(43)

$$\left({\displaystyle \frac{\partial }{\partial t}}+\Omega \frac{\partial }{\partial \theta }\right)\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \tilde{\psi }}{\partial \xi }\right)-\frac{\partial \tilde{\psi }}{\partial \theta }\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)=0\mathrm{if}\xi \in \left(0,1\right)$$

(44)

$$\frac{\partial \tilde{\psi }}{\partial \xi }=0\mathrm{at}\xi =\pm 1$$

(45)

$${\left(\tilde{\psi }\right)}_{\xi =-0}={\left(\tilde{\psi }\right)}_{\xi =+0},{\left(\frac{1}{N^2}\frac{\partial \tilde{\psi }}{\partial \xi }\right)}_{\xi =-0}={\left(\frac{1}{N^2}\frac{\partial \tilde{\psi }}{\partial \xi }\right)}_{\xi =+0}$$

(46)

where the angular velocity $\Omega$ is given by Equation (41).

We confine ourselves to normal modes, i.e., solutions of the form

$$\tilde{\psi }=\widehat{\psi }(r,\xi ){\mathrm{e}}^{\mathrm{i}k\left(\theta -ct\right)}$$

(47)

where k is the azimuthal wavenumber and c is the angular phase speed. Substitution of Equation (47) into Equations (43)–(46) yields

$$\begin{array}{c}\left(c-\Omega \right)\left[\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \widehat{\psi }}{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \widehat{\psi }}{\partial r}\right)-\frac{k^2}{r^2}\widehat{\psi }\right]\hfill \\ +\widehat{\psi }\left\{\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial \left(r^2\Omega \right)}{\partial r}\right]\right\}=0\mathrm{if}\xi \in \left(-1,0\right)\hfill \end{array}$$

(48)

$$\frac{\partial \widehat{\psi }}{\partial \xi }=0\mathrm{at}\xi =-1$$

(49)

$$\left(c-\Omega \right)\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \widehat{\psi }}{\partial \xi }\right)+\widehat{\psi }\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)=0\mathrm{if}\xi \in \left(0,1\right)$$

(50)

$$\frac{\partial \widehat{\psi }}{\partial \xi }=0\mathrm{at}\xi =1$$

(51)

$${\left(\widehat{\psi }\right)}_{\xi =-0}={\left(\widehat{\psi }\right)}_{\xi =+0},{\left(\frac{1}{N^2}\frac{\partial \widehat{\psi }}{\partial \xi }\right)}_{\xi =-0}={\left(\frac{1}{N^2}\frac{\partial \widehat{\psi }}{\partial \xi }\right)}_{\xi =+0}$$

(52)

Together with the requirements that $\widehat{\psi }$ is analytic at $r=0$ and decays as $r\to \infty$, Equations (48)–(52) form an eigenvalue problem, where $\widehat{\psi }(r,\xi )$ is the eigenfunction and c is the eigenvalue. If $\mathrm{Im}c>0$, the vortex characterized by $\Omega (r,\xi )$ and $N\left(\xi \right)$ is unstable.

It can be verified by direct substitution that the solution of the upper-layer problem of Equations (50) and (51) is

$$\widehat{\psi }=a\left(r\right)\left(c-\Omega \right)\mathrm{if}\xi \in \left(0,1\right)$$

(53)

where $a\left(r\right)$ is an arbitrary function. Substituting Equation (53) into the matching conditions of Equation (52) and eliminating $a\left(r\right)$, one obtains

$${\left(\frac{1}{N^2}\frac{1}{\psi }\frac{\partial \widehat{\psi }}{\partial \xi }\right)}_{\xi =-0}=-{\left(\frac{1}{N^2}\frac{1}{c-\Omega }\frac{\partial \Omega }{\partial \xi }\right)}_{\xi =+0}$$

(54)

Due to the inter-layer matching conditions of Equations (40) and (41), “$\xi =+0$” on the right-hand side of Equation (54) can be replaced with “$\xi =-0$”. Replacing “$-0$” with “0”on both sides, one then obtains

$$\frac{\partial \widehat{\psi }}{\partial \xi }=-\frac{1}{c-\Omega }\frac{\partial \Omega }{\partial \xi }\psi \mathrm{at}\xi =0$$

(55)

Thus, the original two-layer eigenvalue problem has been reduced to just the lower-layer problem, comprising Equation (48) and the boundary conditions of Equations (49) and (55).

Finally, we observe that the condition of Equation (55) parametrizes the effect of the upper layer on the stability of the vortex—which does not imply, however, that the upper layer is driven by the lower layer or some such. Equation (55) only applies to a special class of solutions (linear normal modes), whereas the full upper-layer Equation (36) describes independent dynamics and is influenced by the lower layer only through the matching conditions of Equation (52).

4.2. The Equivalent of Rayleigh’s Theorem

Since Lord Rayleigh’s 1880 paper [25], numerous theorems have been proved for parallel and circular flows, relating their stability properties to the across-the-streamlines profile of the (potential) vorticity field. In application to QG parallel flows (a case similar to ours), such a theorem has been proved in [26].

To derive an equivalent of Rayleigh’s Theorem for the problem at hand, consider

$${\int }_0^{\infty }{\int }_{-1}^0(48)\times \frac{r{\widehat{\psi }}^{*}}{c-\Omega }\mathrm{d}\xi \mathrm{d}r,$$

which yields

$$\begin{array}{c}{\int }_0^{\infty }{\int }_{-1}^{0}r{\widehat{\psi }}^{*}\left[\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \widehat{\psi }}{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial \widehat{\psi }}{\partial r}\right)-\frac{k^2}{r^2}\widehat{\psi }\right]\mathrm{d}\xi \mathrm{d}r\hfill \\ +{\int }_0^{\infty }{\int }_{-1}^0\frac{r{\left|\widehat{\psi }\right|}^2}{c-\Omega }\left\{\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial \left(r^2\Omega \right)}{\partial r}\right]\right\}\mathrm{d}\xi \mathrm{d}r=0\hfill \end{array}$$

Integrating by parts the first and second terms of the first integral with respect to $\xi$ and r, respectively, taking into account (49) and (55), and taking imaginary part of the resulting equality, we obtain

$$\left(\mathrm{Im}c\right){\int }_0^{\infty }\left[{\left(\frac{r{\left|\widehat{\psi }\right|}^2}{{\left|c-\Omega \right|}^2}\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)}_{\xi =0}-{\int }_{-1}^0\frac{r{\left|\widehat{\psi }\right|}^2}{{\left|c-\Omega \right|}^2}Q\mathrm{d}\xi \right]\mathrm{d}r=0$$

(56)

where

$$Q=\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial \left(r^2\Omega \right)}{\partial r}\right]$$

is, physically, the radial gradient of the vortex’s PV divided by r.

Identity (56) implies the following stability criterion:

Let Q and ${\left(\partial \Omega /\partial \xi \right)}_{\xi =0}$ be sign-definite and of opposite signs for all $r\in \left(0,\infty \right)$ and $\xi \in \left(-1,0\right)$. Then (56) implies that either $\mathrm{Im}c=0$ or the stability eigenvalue problem does not admit any solutions. Either way, the corresponding vortex is stable.

4.3. Examples of Stable Vortices

We consider sign-definite opposite-sign continuous functions ${\Omega }_0^{\prime }\left(r\right)$ and $Q_0(r,\xi )$, such that

$${\int }_0^{\infty }r{\Omega }_0^{\prime }\left(r\right)\mathrm{d}r<\infty ,{\int }_0^{\infty }{\int }_{-1}^{0}r{Q}_0(r,\xi )\mathrm{d}\xi \mathrm{d}r<\infty$$

Next, we consider the following boundary-value problem for an unknown function $\Omega (r,\xi )$:

$$\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial \left(r^2\Omega \right)}{\partial r}\right]=Q_0(r,\xi )$$

(57)

$$\frac{\partial \Omega }{\partial \xi }=0\mathrm{at}\xi =-1$$

(58)

$$\frac{\partial \Omega }{\partial \xi }={\Omega }_0^{\prime }\left(r\right)\mathrm{at}\xi =0$$

(59)

Together with the requirements that $\Omega (r,\xi )$ is regular at $r=0$ and decays as $r\to \infty$, the boundary-value problem of Equations (57)–(59) has a unique solution describing the lower-layer part of a stable vortex.

Its upper-layer part should be smooth for $\xi >0$ and should satisfy the matching conditions at the interface,

$${\left(\Omega \right)}_{\xi =+0}={\left(\Omega \right)}_{\xi =-0},{\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)}_{\xi =+0}={\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)}_{\xi =-0},$$

but can be arbitrary otherwise.

4.4. Discussion

(1) To understand the difference between the standard PV-based stability criterion and that proved above, we imagine that the upper-layer Equation (50) is treated the same way as the lower-layer Equation (48), i.e., instead of solving Equation (50), we multiply it by $r{\widehat{\psi }}^{*}/\left(c-\Omega \right)$ and integrate, etc. The resulting criterion (to be referred to as C2) states that a vortex is stable if

$$\begin{array}{cc}\hfill Q=\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)+\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{1}{r}\frac{\partial \left(r^2\Omega \right)}{\partial r}\right]& \mathrm{does}\mathrm{not}\mathrm{change}\mathrm{sign}\mathrm{if}\xi \in \left(-1,0\right),\phantom{{\displaystyle \underset{a}{\int }}}\hfill \\ \hfill Q=\frac{\partial }{\partial \xi }\left(\frac{1}{N^2}\frac{\partial \Omega }{\partial \xi }\right)& \mathrm{does}\mathrm{not}\mathrm{change}\mathrm{sign}\mathrm{if}\xi \in \left(0,1\right),\hfill \end{array}$$

and $Q_{\xi >0}$ and $Q_{\xi <0}$ are of the same sign.

Thus, C2 guarantees stability only if the vortex’s PV gradient is sign-definite in the whole bulk of the fluid (that Q has a "reduced" form in the upper layer is unimportant). Clearly, this is a much smaller class of vortices than that described by our criterion, as the latter does not restrict the $\xi >0$ part of the vortex.

We have obtained a wider class of stable vortices due to the explicit solution of Equation (53), which would be impossible to derive in, for example, the standard QG formulation.

(2) The most unexpected feature of our criterion is that it allows the $\xi >0$ part of a stable vortex to be arbitrary.

This result—no matter how surprising—agrees with that obtained in [14,15] for the two-layer model. According to these papers, if the lower-layer PV is constant (in terms of the present work, this corresponds to $Q_0=0$), all vortices with $B{u}_1\ll 1$ are stable, regardless of all other upper-layer characteristics. This result has been obtained asymptotically for QG vortices with a thin upper layer, and confirmed numerically for the exact models: both the QG and the primitive.

One concludes that a thick layer with a stable flow in it can "absorb" any potential instability of the upper-layer flow.

5. Concluding Remarks

Summarizing the present work, we state that its main result is the wide class of stable vortices presented in Section 4.3. It generalizes the class discovered in [14,15] and has the potential to explain the observed longevity of oceanic eddies.

It remains to list the caveats and disclaimers that our conclusions are subject to.

(1)

The applicability of our asymptotic results is restricted by the use of the same velocity scale for both active and passive layers, whereas a model targeting oceanic eddies should allow the flow in the former to exceed that in the latter by an order of magnitude.

Such a modification (scaling down the lower-layer velocity) should not change our main conclusion, however. As shown for the two-layer approximation of stratification [12,13,14,15], even a weak flow in the passive layer is capable of stabilizing the vortex, and the same is suggested by preliminary estimates for continuous stratification (this work is currently in progress).

(2)

All of our results are based on the linear stability of vortices, and, thus, neglect nonlinearity. In principle, nonlinear effects can provide an alternative mechanism of vortex stabilization: unstable perturbations can saturate at a certain amplitude [3,4,27,28], or unstable (circular) vortices can perhaps evolve into stable elliptic vortices [29] or stable tripoles [30]. So far, however, these models have not been tested for the parameters of real oceanic eddies.

(3)

Our results apply to vortices for which ageostrophic effects are weak, and so they have been neglected—if they were taken into account, they might potentially cause a weak, higher-order instability. We note that this does not occur for two-layer vortices, for which the QG model [9,14] and the primitive equations [10,15] yield remarkably similar results. Nevertheless, to be on the safe side, the stable vortices found in this paper should be re-examined using the primitive equations.

However, even if a weak ageostrophic instability is found for the vortices presented in this paper, the longevity of oceanic eddies can still be accounted for if the instability’s e-folding time (the reciprocal of the growth rate) is of the order of the eddies’ characteristic lifetime.