Vortex Stability in the Thermal Quasi-Geostrophic Dynamics

1. Introduction

Vortices are an essential feature of planetary fluid flows. Their stability has often been studied using the rotating shallow water (RSW) model with one or several homogeneous layers (i.e. with a mean vertical stratification) [1,2,3,4,5]. In particular, these studies showed that vortices with horizontally or vertically sheared velocity profiles can be sensitive to barotropic or baroclinic instabilities. In their nonlinear evolution, strongly unstable vortices break into several fragments. Only for weakly unstable vortices does the perturbation stabilize non-linearly, leading to a vortex multipole. Vortex stability has also been studied in the quasi-geostrophic model, a subset of the RSW equations, valid for low frequency dynamics [6,7,8,9]. These studies also showed that very unstable vortices would break irreversibly into dipoles, while weakly unstable vortices could rearrange as rotating multipoles.

However, neither the RSW model nor the usual quasi-geostrophic model, are appropriate in the presence of a mean horizontal stratification. Therefore, an extension of RSW models to horizontally inhomogeneous layers was derived by [10]. This model was called ILPE (Inhomogeneous Layer Primitive Equation model) — or TRSW (Thermal Rotating Shallow-Water model) for the specific case of the $I{L}^0$ model — in subsequent papers. In his 1993 paper, P. Ripa showed that the pseudo-energy or the pseudo-momentum of flow could not be ascertained to be sign-definite for general flows and therefore that no general criterion for formal stability could be obtained in the ILPE model. This in no way restricts the usefulness of the ILPE nor of the TRSW model. Indeed, they were used to model the upper ocean dynamics [11,12,13] and they were also developed to include moisture for application to the atmosphere [14,15].

The TRSW model equations being difficult to handle, several authors derived a low-Rossby number (low-frequency) version of these equations, called TQG (Thermal Quasi-Geostrophy). Firstly, P. Ripa derived a low-frequency approximation of the ILPE model that he called ILQG [16]. The equations of this model are those later called thermal quasi-geostrophic (TQG) equations, also derived and studied by [17] and by [18]. A numerical analysis of the TQG equations was achieved by [19]. [20] derived a criterion for the nonlinear saturation of thermal instabilities in this model.

[21] studied numerically the stability of jets and vortices in the TRSW model and showed that small-scale features appear as the instability develops. The complexity of the TRSW model renders difficult any analytical work on jet or vortex stability. Therefore, we use the TQG model here, to analyze the instability of circular vortices. After recalling the model equations (Section 2), we provide analytical solutions for linear stability in a vortex core, and then a general stability criterion for circular flows (Section 3). We also compute the growth rates and unstable waves for several velocity and buoyancy distributions numerically. In Section 4, we use a numerical model to classify the nonlinear regimes of thermally unstable vortices in the TQG dynamics. We analyze the finite-time evolution of these vortices with a Fourier analysis and the frontogenesis function. Finally, we conclude this study.

2. Model

2.1. Equations

We mostly follow the notations introduced by [19] to express the two equations of the Thermal Quasi-Geostrophic model as:

$$\begin{array}{ccc}\hfill {\partial }_t\omega +\mathbf{u}·\nabla \omega & =& \mathbf{u}·\nabla b-{\mathbf{u}}_{{\mathbf{h}}_{\mathbf{1}}}·\nabla b,\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill {\partial }_{t}b+\mathbf{u}·\nabla b& =& 0,\hfill \end{array}$$

(2)

where

$$\omega =\left({\nabla }^2-{\displaystyle \frac{1}{R_d^2}}\right)\psi +f,$$

(3)

is the potential vorticity, $\psi$ is the streamfunction, $R_d=\sqrt{gH}/f$ is the deformation radius (taken as unity by [19]). In these notations, g is gravity, H is the maximal thickness of the fluid layer at rest, and f is the Coriolis parameter, chosen constant here. In these equations, b is the buoyancy (which depends directly on temperature via the equation of state of the fluid, assuming no other active tracer). All variables $\omega ,\psi ,b$ depend only on the horizontal coordinates, being vertically averaged.

In the equations above, the velocities are given by the following relations

$$\begin{array}{ccc}\hfill \mathbf{u}& =& \mathbf{k}\times \nabla \psi ,\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill {\mathbf{u}}_{{\mathbf{h}}_{\mathbf{1}}}& =& \mathbf{k}\times \nabla h_1,\hfill \end{array}$$

(5)

where $h_1$ is the bottom height above the maximal depth and $\mathbf{k}$ the vertical unit vector.

2.2. Flow and buoyancy distributions

As we study vortices here, the mean velocity and buoyancy distributions will be chosen circular, and thus dependent on the radius r only, $\overline{\psi }\left(r\right),\overline{b}\left(r\right)$.

In particular, to derive simple analytical solutions in SubSection 3.1, we will use

$$\overline{\psi }\left(r\right)={\Omega }_0{\displaystyle \frac{r^2}{2}},\overline{b}\left(r\right)=B_0{\displaystyle \frac{r^2}{2}},$$

which is an acceptable form for the core of oceanic vortices, see for instance [22,23]. We define $\overline{V}=d\overline{\psi }/dr$ the tangential velocity of the mean flow (its radial component is null), and $\overline{\Omega }=\overline{V}/r$ the corresponding mean rotation rate. Bottom topography will be included only locally in this study.

For the derivation of the linear instability criterion (in SubSection 3.2 and Section 3.3), a general form of these functions will be kept. Finally, for the numerical study of both linear and nonlinear evolutions of unstable vortices (SubSection 3.4 and Section 4), Gaussian or power-exponential profiles of mean streamfunction and buoyancy will be considered.

Power exponential vortices have the following profile of vorticity and of buoyancy

$$\begin{array}{ccc}\hfill \overline{\omega }\left(r\right)& =& {\omega }_0\left[1-{\displaystyle \frac{\alpha }{2}}{(r/r_0)}^{\alpha }\right]exp(-{(r/r_0)}^{\alpha }),\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill \overline{B}\left(r\right)& =& B_{0}exp(-{(r/r_0)}^{\alpha }),\hfill \end{array}$$

(7)

with ${\omega }_0=2V_0/r_0$, $V_0$ being the vortex azimuthal velocity and $r_0$ its radius. $B_0$ is the vortex buoyancy.

For the linear instability analysis, angular normal modes will be used

$${\psi }^{\prime }(r,\theta ,t)=\varphi \left(r\right)exp\left[il(\theta -ct)\right],b^{\prime }(r,\theta ,t)=b\left(r\right)exp\left[il(\theta -ct)\right],$$

where l is the angular mode (or azimuthal mode), and the complex value of $c=c_r+i{c}_i$ contains the phase speed of the perturbation (the real part of c) and the growth rate of the perturbation ($l{c}_i$). With such perturbations, the linearized equations of vorticity and of buoyancy are

$$\begin{array}{ccc}\hfill (\overline{\Omega }-c){\omega }^{\prime }-\overline{\Omega }{b}^{\prime }-{\displaystyle \frac{{\psi }^{\prime }}{r}}{\displaystyle \frac{d}{dr}}(\overline{\omega }-\overline{b})& =& 0,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill (\overline{\Omega }-c)b^{\prime }-{\displaystyle \frac{{\psi }^{\prime }}{r}}{\displaystyle \frac{d\overline{b}}{dr}}& =& 0,\hfill \end{array}$$

(9)

in the absence of bottom topography ($h_1=0$).

2.3. Linear and nonlinear numerical models

In SubSection 3.4, when using normal mode perturbations and when discretizing the radius over N grid points ($r_1...r_N$), the linear equations 8, become a $2N\times 2N$ matrix system, in $\varphi$ and b. This system is solved by diagonalization, c being the eigenvalue and $X(\varphi \left(r_1\right)...\varphi \left(r_N\right),b\left(r_1\right)...b\left(r_N\right))$ being the corresponding eigenvector. We thus obtain $AX=cBX$, where the A and B matrices depend on $\overline{\Omega },\overline{\omega },\overline{b}$.

In Section 4, the numerical simulations with the system of nonlinear equations 1, will be carried out with a pseudo-spectral code of these equations, on a bi-periodic horizontal grid. The number of grid points is 512 horizontally. The code uses a mixed Euler-leapfrog scheme in time which is conservative in vorticity and buoyancy, and which does not separate odd and even solutions in time-steps. The horizontal derivatives are performed in spectral space, using fast Fourier transforms to recover physical variables. Very weak bi-harmonic viscosity is added to equations 1, which does not alter the physical results (tests have been performed with different small values of viscosity).

3. Linear stability

3.1. A simple analytical solution

Firstly, we look for an analytical solution to equations 8, ; this is possible only for simple mean streamfunction and buoyancy profiles. Here we choose such profiles, which are realistic for the core of vortices. It is not possible to solve analytically the problem for the whole vortex, but in the core resides most of the mean kinetic (or potential) energy which feeds the perturbation.

Choosing

$$\overline{\psi }\left(r\right)={\Omega }_0{\displaystyle \frac{r^2}{2}},\overline{b}\left(r\right)=B_0{\displaystyle \frac{r^2}{2}},\overline{\Omega }={\Omega }_0,{\displaystyle \frac{1}{r}}{\displaystyle \frac{d\overline{b}}{dr}}=B_0$$

and inserting these expressions into equations 8, (with flat bottom dynamics), leads to

$$\begin{array}{ccc}\hfill ({\Omega }_0-c){\omega }^{\prime }-{\Omega }_0{b}^{\prime }-{\psi }^{\prime }\left({\displaystyle \frac{1}{r}}{\displaystyle \frac{d\overline{\omega }}{dr}}-B_0\right)& =& 0,\hfill \end{array}$$

(10)

$$\begin{array}{ccc}\hfill ({\Omega }_0-c)b^{\prime }-B_0{\psi }^{\prime }& =& 0.\hfill \end{array}$$

(11)

The mean potential vorticity is $\overline{\omega }={\Omega }_0(2-r^2/\left(2R_d^2\right))$ and therefore

$${\displaystyle \frac{1}{r}}{\displaystyle \frac{d\overline{\omega }}{dr}}=-{\displaystyle \frac{{\Omega }_0}{R_d^2}}.$$

(12)

With this, the linear equations are

$$\begin{array}{ccc}\hfill ({\Omega }_0-c){\omega }^{\prime }-{\Omega }_0{b}^{\prime }+{\psi }^{\prime }\left({\displaystyle \frac{{\Omega }_0}{R_d^2}}+B_0\right)& =& 0,\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill ({\Omega }_0-c)b^{\prime }-B_0{\psi }^{\prime }& =& 0.\hfill \end{array}$$

(14)

The last step relates ${\omega }^{\prime }$ and ${\psi }^{\prime }$. This is achieved by considering eigen-functions of the Laplacian operator in polar coordinates

$${\psi }^{\prime }(r,\theta ,t)={\psi }_0{J}_l(r/R_d)exp\left[il(\theta -ct)\right],$$

(15)

which leads to ${\omega }^{\prime }=-2{\psi }^{\prime }/R_d^2$ (for simplicity, we have not specified here that we take the real part of the expression above; therefore ${\psi }_0$ is complex).

Similarly, we set

$$b^{\prime }(r,\theta ,t)=b_0{J}_l(r/R_d)exp\left[il(\theta -ct)\right].$$

(16)

This provides a $2\times 2$ homogeneous system in ${\psi }_0,b_0$ whose determinant must vanish. This determinant is

$$D=(2c-{\Omega }_0)({\Omega }_0-c)-c{B}_0{R}_d^2=0.$$

(17)

Instability will occur in such a vortex core if

$$\Delta ={({\Omega }_0-B_0{R}_d^2)}^2-4{\Omega }_0{B}_0{R}_d^2<0.$$

(18)

Thus, a necessary condition for instability in this case is ${\Omega }_0{B}_0>0$.

We can analyze condition $\Delta <0$ more fully. Setting $X={\Omega }_0,Y=B_0{R}_d^2$, this condition can also be written as $(X-Y+2\sqrt{XY})(X-Y-2\sqrt{XY})<0$. Some short algebra indicates that this is possible only if $X-Y<2\sqrt{XY}$ and $X-Y+2\sqrt{XY}>0$.

When adding topography to the problem, we can choose a parabolic seamount ${\overline{h}}_1={\Omega }_1{r}^2/2$, so that ${\overline{V}}_1/r={\Omega }_1$. Then we have a new condition for instability

$$\Delta ={({\Omega }_0-B_0{R}_d^2)}^2-4({\Omega }_0-{\Omega }_1)B_0{R}_d^2<0.$$

(19)

Since ${\Omega }_1>0$, a seamount stabilizes the flow if $B_0>0$, and conversely.

3.2. Stability criterion

After this first simple solution, we look for a general stability criterion for TQG vortices. Now we do not choose a specific form of $\overline{\psi }$, nor of $\overline{b}$ any more. Equation 8 could be amenable to a Rayleigh-type criterion by multiplying it by ${\psi }^{\prime *}$ (the complex conjugate of ${\psi }^{\prime }$), but then the coupled term in ${\psi }^{\prime *}{b}^{\prime }$ should be eliminated; this implies to multiply equation by the proper quantity for a further substitution. This is done in SubSection 3.3. Here we choose to express $b^{\prime }$ in terms of ${\psi }^{\prime }$ using equation , without multiplying equations 8, by complex conjugates of ${\psi }^{\prime }$ or $b^{\prime }$ immediately.

Setting

$${\psi }^{\prime }=\varphi \left(r\right)exp\left(il(\theta -ct)\right),b^{\prime }=b\left(r\right)exp\left(il(\theta -ct)\right)$$

and using the expression of b with respect to $\varphi$ from equation , equation 8 becomes

$$(\overline{\Omega }-c)\left[\ddot{\varphi }+{\displaystyle \frac{\dot{\varphi }}{r}}-{\displaystyle \frac{l^2}{r^2}}\varphi -{\displaystyle \frac{\varphi }{R_d^2}}\right]-\varphi {\displaystyle \frac{\overline{\Omega }\dot{\overline{b}}}{r(\overline{\Omega }-c)}}-\varphi {\displaystyle \frac{[\dot{\overline{\omega }}-\dot{\overline{b}}]}{r}}=0,$$

(20)

where $\dot{f}\left(r\right)=df/dr$.

This equation resembles the Taylor-Goldstein equation. This can be understood since they both result from a vorticity and a buoyancy equation (in different planes, horizontal or vertical). Thus, it is appropriate to use the change of variable $\varphi \left(r\right)=\sqrt{\overline{\Omega }-c}\chi \left(r\right)$. Setting $\alpha =\sqrt{\overline{\Omega }-c}$, some algebra leads to the equation

$${\displaystyle \frac{1}{r}}{\displaystyle \frac{d}{dr}}\left[r(\overline{\Omega }-c)\dot{\chi }\right]-{\displaystyle \frac{l^2}{r^2}}{\alpha }^2\chi -{\displaystyle \frac{1}{R_d^2}}{\alpha }^2\chi +\chi {\displaystyle \frac{\ddot{\overline{\Omega }}}{2}}-{\displaystyle \frac{\chi }{r}}(\dot{\overline{\omega }}-\dot{\overline{b}})+{\displaystyle \frac{\dot{\overline{\Omega }}\chi }{2r}}-{\displaystyle \frac{\chi {\left(\dot{\overline{\Omega }}\right)}^2}{4{\alpha }^2}}-{\displaystyle \frac{\chi }{r}}{\displaystyle \frac{\overline{\Omega }\dot{\overline{b}}}{{\alpha }^2}}=0,$$

(21)

where $\dot{\chi }$ is the derivative of $\chi$.

Only now do we multiply this expression by ${\chi }^{*}$ which is the complex conjugate of $\chi$ and we integrate the equation over r as $\int rdr$. This leads to

$$-\int (\overline{\Omega }-c)\left[\right|\dot{\chi }{|}^2+\left({\displaystyle \frac{l^2}{r^2}}+{\displaystyle \frac{1}{R_d^2}}\right){\left|\chi \right|}^2{]rdr+\int |\chi |}^2\left[{\displaystyle \frac{\ddot{\overline{\Omega }}}{2}}+{\displaystyle \frac{\dot{\overline{\Omega }}}{2r}}-{\displaystyle \frac{1}{r}}(\dot{\overline{\omega }}-\dot{\overline{b}})\right]rdr-$$

$$\int {\displaystyle \frac{\frac{{\left(\dot{\overline{\Omega }}\right)}^2}{4}+\frac{\overline{\Omega }\dot{\overline{b}}}{r}}{\overline{\Omega }-c}}{\left|\chi \right|}^{2}rdr=0.$$

(22)

We can take the imaginary part of this expression; the central term will vanish and we will be left with

$$c_i\left[\int \left[\right|\dot{\chi }{|}^2+\left({\displaystyle \frac{l^2}{r^2}}+{\displaystyle \frac{1}{R_d^2}}\right)\left|\chi {|}^2\right]rdr\right]-c_i\left[\int {\displaystyle \frac{\frac{{\left(\dot{\overline{\Omega }}\right)}^2}{4}+\frac{\overline{\Omega }\dot{\overline{b}}}{r}}{|\overline{\Omega }{-c|}^2}}{\left|\chi \right|}^{2}rdr\right]=0.$$

(23)

For $c_i$ to be non-zero, the second integral must be positive. Therefore, we must have

$$C\left(r\right)={\displaystyle \frac{-\overline{\Omega }d\overline{b}/dr}{{\left(\dot{\overline{\Omega }}\right)}^{2}r}}<1/4,$$

(24)

at some r, for instability; this is a necessary condition.

We can also write this inequality as

$$\overline{\Omega }d\overline{b}/dr>-r({\dot{\overline{\Omega })}}^2/4.$$

(25)

Applying this inequality to the case of SubSection 3.1, we obtain ${\Omega }_0{B}_0>0$, for instability, which is the condition previously found.

We can generalize this criterion to the presence of bottom topography $h_1\left(r\right)$. Defining ${\overline{\Omega }}_1\left(r\right)=(1/2r)d{\overline{h}}_1/dr$, we need to have

$$C_1\left(r\right)={\displaystyle \frac{[{\overline{\Omega }}_1-\overline{\Omega }]d\overline{b}/dr}{{\dot{\left(\overline{\Omega }\right)}}^{2}r}}<1/4,$$

(26)

at some r, for instability. Again, we can see that a seamount is stabilizing for positive $d\overline{b}/dr$.

3.3. A second instability criterion

As mentioned in the previous subsection, we can multiply equation 8 by ${\psi }^{\prime *}$, the complex conjugate of ${\psi }^{\prime }$ and divide it by $(\overline{\Omega }-c)$, leading to

$${\omega }^{\prime }{\psi }^{\prime *}-{\displaystyle \frac{\overline{\Omega }}{(\overline{\Omega }-c)}}{b}^{\prime }{\psi }^{\prime *}-{\displaystyle \frac{|{\psi }^{\prime }{|}^2}{r}}{\displaystyle \frac{1}{(\overline{\Omega }-c)}}{\displaystyle \frac{d}{dr}}[\overline{\omega }-\overline{b}]=0.$$

(27)

we also multiply the complex conjugate of equation by $b^{\prime }$ and divide it by $(\overline{\Omega }-c^{*})$, leading to

$$|b^{\prime }{|}^2-{\displaystyle \frac{d\overline{b}/dr}{(\overline{\Omega }-c^{*})}}{b}^{\prime }{\psi }^{\prime *}=0.$$

(28)

Finally, we eliminate $b^{\prime }{\psi }^{\prime *}$ between these two equations, and integrate $\int rdr$ to obtain

$$K^{\prime }+P^{\prime }=-\int {\displaystyle \frac{rdr}{(\overline{\Omega }-c)}}[{\displaystyle \frac{\overline{\Omega }(\overline{\Omega }-c^{*})}{d\overline{b}/dr}}|b^{\prime }{|}^2+{\displaystyle \frac{|{\psi }^{\prime }{|}^2}{r}}{\displaystyle \frac{d}{dr}}(\overline{\omega }-\overline{b})],$$

(29)

where

$$(K^{\prime },P^{\prime })=\int rdr\left(\right|{\partial }_r{\psi }^{\prime }{|}^2+l^2|{\psi }^{\prime }{|}^2/r^2,|{\psi }^{\prime }{|}^2/R_d^2),$$

are (twice) the kinetic, potential energies of the perturbation.

Let us now take the imaginary part of this equation; it is

$$c_i\int rdr[{\displaystyle \frac{1}{|\overline{\Omega }{-c|}^2}}{\displaystyle \frac{2\overline{\Omega }(\overline{\Omega }-c_r)}{d\overline{b}/dr}}|b^{\prime }{|}^2+{\displaystyle \frac{|{\psi }^{\prime }{|}^2}{r}}{\displaystyle \frac{d}{dr}}(\overline{\omega }-\overline{b})]=0.$$

(30)

For the vortex to be unstable, $c_i$ must be non-zero (positive), and therefore the two terms in the integral must be opposite-signed. Thus, a necessary condition for instability is

$$2\overline{\Omega }(\overline{\Omega }-c_r){\displaystyle \frac{d\overline{b}}{dr}}\mathrm{and}{\displaystyle \frac{d}{dr}}[\overline{\omega }-\overline{b}]\mathrm{must}\mathrm{be}\mathrm{opposite}\mathrm{signed}.$$

(31)

Since this criterion involves $c_r$, it is less straightforward to apply.

We have not found how to substitute $c_r$ by some $\overline{\Omega }\left(r_c\right)$ rigorously, as in the Fjortoft criterion. Obtaining $c_r$ exactly requires solving the linear instability problem (numerically), which is not possible for all mean flow and buoyancy distributions. A rough estimate of $c_r$ could be computed from our analytical solution at marginality.

3.4. Numerical analysis of linear stability

Now we choose a realistic mean distribution of $\overline{V},\overline{b}$ to discretize equations 8, over the N steps in r and to solve the linear instability problem numerically. We start with Gaussian mean distributions of streamfunction and of buoyancy

$$\begin{array}{ccc}\hfill \overline{V}\left(r\right)& =& V_0(r/r_0)exp(-{(r/r_0)}^2),\hfill \end{array}$$

(32)

$$\begin{array}{ccc}\hfill \overline{b}\left(r\right)& =& B_{0}exp(-{(r/r_0)}^2).\hfill \end{array}$$

(33)

For simplicity, we set $V_0=1,r_0=1$. This provides a length and a time scale. We also set $F_1=1/R_d=0$ unless otherwise stated.

As a reference, we solve equations 8, with $B_0=0,b^{\prime }=0$; this is the 2D incompressible flow case that is well known [24,25]). For the TQG case, we set $B_0=0.1$.

Figure 1a shows the growth rate $\sigma =l{c}_i$ versus the azimuthal wavenumber l for both cases. In the 2D flow, only $l=2$ is unstable for Gaussian vortices. In the TQG flow, the Gaussian vortex is linearly unstable for all wavenumbers from $l=1$ to $l=5$. In TQG, there is no universally stable solution for $l=1$, contrary to 2D [24]. Furthermore, all growth rates are much larger in TQG than in 2D (twice as large for $l=2$). From this, we can conclude that TQG vortices are not only more unstable than their 2D counterparts, but also that they will develop smaller-scale perturbations. The origin of these differences between the two models can be searched in equations 8, . Firstly, it is clear that if we set $d\overline{b}/dr=0$ in both equations, then $b^{\prime }=0$, and equation 8 becomes the 2D linear equation. Therefore, a non-zero mean buoyancy gradient is necessary to account for this difference. It is also of interest to see that if we cancel the mean buoyancy gradient in the second equation only, the growth rates in this case are comparable to those of 2D as long as $B_0/V_0\ll 1$. In this case, like-signed $B_0$ and $V_0$ increase the mean vorticity/buoyancy gradient. More interesting is the cancellation of the $\overline{\Omega }{b}^{\prime }$ term only, in equation 8; this renders the first equation 2D-like (independent from $b^{\prime }$); equation then only provides $b^{\prime }$ as a function of ${\psi }^{\prime }$ and the growth rates are those of 2D, again in the limit $B_0/V_0\ll 1$. Therefore, clearly, the term $\overline{\Omega }{b}^{\prime }$ is responsible for the difference between the two dynamics, TQG on the one hand, and 2D with a passive b tracer, on the other hand.

For the same (Gaussian) flow, we study now the influence of the ratio $B_0/V_0$ on the growth rate $\sigma =l{c}_i$ for $l=2$ (see Figure 1b). Clearly, in the TQG model, the coupling of the vorticity equation with the buoyancy equation destabilizes the flow for small values of this ratio (compared to the 2D model). Increasing this ratio to unity then decreases the growth rates. A maximal growth rate is obtained for $B_0/V_0\sim 0.3$.

Again for the Gaussian vortex, with $l=2$ and with $B_0=0.1,V_0=1.0$, we test the sensitivity of growth rates to free surface effects. We call $F_1=1/R{d}^2$ and plot on Figure 1c, the growth rates with respect to $F_1$ both in 2D and in TQG. In TQG, the growth rates of horizontal shear instability decrease when $F_1$ increases, as in 2D flows [26].

Finally, still for the Gaussian vortex, with $l=2$ and with $B_0=0.1,V_0=1.0$, we test the sensitivity of growth rates to bottom topography. This bottom topography is chosen via ${\Omega }_1\left(r\right)={\Omega }_{1}exp(-r^2)$ (see its formal expression in SubSection 3.2). Figure 1d confirms that a seamount stabilizes the flow if $B_0>0$.

In 2D flows, it is well known that the growth rates of the perturbation increase with the mean vorticity shear of the vortex. To this aim, we now use the following family of profiles, also known as $\alpha$-exponential profiles [27], for the velocity

$$\begin{array}{ccc}\hfill \overline{V}\left(r\right)& =& V_0(r/r_0)exp(-{(r/r_0)}^{\alpha }),\hfill \end{array}$$

(34)

$$\begin{array}{ccc}\hfill \overline{b}\left(r\right)& =& B_{0}exp(-{(r/r_0)}^{\alpha }),\hfill \end{array}$$

(35)

Figure 2a shows the growth rates versus $\alpha$ in the 2D case ($B_0=b^{\prime }=0$), and Figure 2b in the TQG case.

As usual for 2D dynamics, we observe that $\sigma$ increases with $\alpha$ for $l=2,3$; $\alpha$ is not large enough in the chosen range to allow the growth of modes $l=4,5$. We recover the critical value of $\alpha \sim 1.9$ for $l=2$ as mentioned in [27]. Barotropic instability occurs only for steep enough vortices, all the more so as the waves are short.

On the contrary, in TQG dynamics, vortices with "flat" vorticity profiles are unstable to both long and short waves. A small growth rate $\sigma \sim 0.03$ exists for $l=1-5$ over a wide range of values of $\alpha$. The $\sigma$ curves for $l=2,3$ emerge from this minimal growth rate; it must be noted that for $\alpha$ large enough ($\sim 3-4$), the growth rates in the 2D and TQG are comparable for $l=2,3$.

Physically, one must exert caution for $\alpha \le 1.0$, though the results are mathematically correct. For $\alpha =1.0$, we have $\overline{\omega }=2V_0/r_0[2-r/r_0]exp(-r/r_0)$. Vorticity behaves as $exp(-r)$ near the center and, in a two-dimensional representation, exhibits a cusp at the center. Therefore, it is recommended to consider only cases with $\alpha >1$.

4. Nonlinear evolution

After the linear analysis of stability, we study numerically the nonlinear regimes of perturbed unstable vortices. Firstly, we study a reference case with Gaussian distributions of buoyancy and streamfunction, with a given ratio between their amplitude, and an elliptical perturbation. Then we will vary all these parameters, and briefly study the influence of a finite deformation radius and of circular topography.

4.1. Reference case

Our reference case has

$$\psi ={\psi }_{0}exp(-r^2/r_0^2),b=b_{0}exp(-r^2/r_0^2),b_0/{\psi }_0=0.3$$

with only an elliptical perturbation initially (Figure 3a-b at time $t_0=0$). With these conditions, we ran the TQG model and the decoupled 2D model in parallel. Figure 3 shows the time evolution of vorticity and of buoyancy in the two cases. Initially, the decoupled 2D model dynamics creates filaments during the reorganization of the vortex periphery, yet fewer than the TQG model (see fig.3a-b at time $t_1=100$). At finite times, the TQG model creates many small scale features than the decoupled model, in particular radially. This is not the case of the 2D decoupled model (same figures at time $t_2=200$).

To interpret this appearance of fronts, we note again that the term $\overline{\Omega }{b}^{\prime }$ plays a key role in the coupling of the vorticity and buoyancy equations. Vorticity and thermal Rossby waves propagate circularly on the radial vorticity and buoyancy gradients. Since these gradients (and the mean flow) are not constant, the buoyancy anomaly of these waves is sheared and forms filaments, Via the coupling, vorticity anomalies are formed; they are also sheared. This is clearly seen in the TQG model vorticity plots (Figure 3a). The buoyancy distribution being more spread out that that of vorticity, filaments are formed mostly around the vortex core. The vorticity anomalies lead to shear flow growth, which increase ${\omega }^{\prime }$ and $b^{\prime }$.

To quantify this process further, we decompose the vorticity anomalies azimuthally in Fourier modes, and radially in Bessel components. This is done via

$$\begin{array}{ccc}\hfill {\omega }^{\prime }(r,\theta ,t)& =& \sum _l{C}_l(r,t)cos\left(l\theta \right)+S_l(r,t)sin\left(l\theta \right),\hfill \end{array}$$

(36)

$$\begin{array}{ccc}\hfill A_l\left(t\right)& =& (2/R^2)\int rdr[C_l^2+S_l^2],\hfill \end{array}$$

(37)

where R is an arbitrary radius (here chosen as $\pi$), large enough to encompass the vortex core and periphery. The modal amplitudes $A_l\left(t\right)$ are plotted on Figure 4a,c. For sake of clarity l was considered continuous on the graph, but only the integer values of l have to be considered. Clearly, mode $l=2$ remains the dominant (and nearly sole) mode in the decoupled evolution, while higher modes are much stronger in the TQG model. They correspond, in physical space, to the small-scale features observed in the vorticity field.

Secondly, for $l=2$, we project $C_l(r,t)$ and $S_l(r,t)$ onto a basis of Bessel functions using their orthogonality.

$$\begin{array}{ccc}\hfill C_{2k}\left(t\right)& =& \int rdr{C}_2(r,t)J_2({\alpha }_{2k}r/R_0),\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill S_{2k}\left(t\right)& =& \int rdr{S}_2(r,t)J_2({\alpha }_{2k}r/R_0),\hfill \end{array}$$

(39)

and

$$A_{2k}\left(t\right)=\sqrt{C_{2k}^2+S_{2k}^2}\left(t\right),$$

(40)

where $R_0$ is the vortex radius (here equal to unity), and ${\alpha }_{2k}$ is the $k^{th}$ root of the Bessel function $J_2$. It is known that for increasing ${\alpha }_{2k}$, the Bessel function $J_2$ oscillates faster radially (corresponding to thinner structures). For sake of clarity k was considered continuous on the graph, but only the integer values of k have to be considered. Figure 4b,d indicate that radially short perturbations (components $k=3,4,5$) grow initially on the vortex in the decoupled 2D dynamics, but they rapidly decay afterwards. On the contrary, short radial perturbations remain of larger amplitude in the TQG dynamics for the whole duration of the simulation.

Now, we compute the two components of the deformation, the normal strain rate ${\sigma }_n={\partial }_{x}u-{\partial }_{y}v$, the shear strain rate ${\sigma }_s={\partial }_{x}v+{\partial }_{y}u$ and, from them, the total strain rate $\sigma =\sqrt{{\sigma }_n^2+{\sigma }_s^2}$. We plot this latter, for $t=100$ in the simulation, on Figure 5a,b for the TQG and the decoupled models. Clearly, the deformation rate is stronger and more spread out in the TQG simulation. Also many more small-scale intricate patterns are seen in this simulation compared with the 2D decoupled model one.

Finally, using the Lagrangian conservation of buoyancy in both TQG and 2D equations, one can derive an equation for the evolution of the gradient of buoyancy (equivalently, of buoyancy)

$${\displaystyle \frac{d\nabla b}{dt}}=-\nabla \mathbf{u}·\nabla b,$$

(41)

which can be multiplied by the adjoint of the buoyancy gradient to yield

$${\displaystyle \frac{1}{2}}{\displaystyle \frac{{d|\nabla b|}^2}{dt}}=-\nabla b·\nabla \mathbf{u}·\nabla b=F.$$

(42)

The frontogenetic tendency, or frontogenesis function, F indicates where the buoyancy (or buoyancy) gradients will amplify or decay, and thus where filaments will grow or disappear. This function is shown for our TQG simulation on Figure 6. The regions shown correspond to the regions of strong buoyancy gradient in Figure 3a at $t=100$. Clearly, buoyancy gradients are formed in the vicinity of the deformed vortex core; they tend to disappear in the periphery.

4.2. Parameter sensitivity

4.2.1. Vorticity profile steepness and relative intensity of the mean buoyancy

Firstly, for an elliptical perturbation ($l=2$), we vary the steepness ($\alpha$) of the mean velocity and buoyancy profiles and the ratio of their intensities ($B_0/V_0$). Here we have

$$\begin{array}{ccc}\hfill \overline{V}\left(r\right)& =& V_0(r/r_0)exp(-{(r/r_0)}^{\alpha }),\hfill \end{array}$$

(43)

$$\begin{array}{ccc}\hfill \overline{b}\left(r\right)& =& B_{0}exp(-{(r/r_0)}^{\alpha }).\hfill \end{array}$$

(44)

Figure 7a indicates that for Gaussian vortices ($\alpha =2$), tripoles are formed by the instability of the circular vortex. Tripole formation is all the less efficient as $B_0/V_0$ grows. For large values of this ratio, the end-state of the nonlinear simulation is a quasi-axisymmetric vortex. Indeed, as mean buoyancy is larger, the vortex is more prone to small-scale instabilities (as was shown by linear instability theory). This renders the satellites of the tripole less coherent. Thus, their robustness is lesser and the shear that they exert on the vortex core is weaker. For large $B_0/V_0$, the satellites are not coherent any more and the tripolar state is only transient.

For steeper vortices, the present of filaments in the tripole satellites and around them becomes more prevalent until the whole vortex periphery is completely turbulent (for growing $\alpha$, with constant $B_0/V_0$; see Figure 7b). Also a mode $l=4$ is observed in the vortex evolution, due both to wave-wave interaction (from the fundamental mode $l=2$) and to the growth of small waves in the TQG model (see also Figure 7b). For steep vortices, the radial gradient of mean vorticity is large enough to lead to vortex breaking into two dipoles (the usual irreversible, fully nonlinear, evolution of strongly unstable vortices with $l=2$ see Figure 7c). Then, nonlinear harmonic formation from $l=2$ is more efficient than the growth of small waves via the TQG dynamics. Again, for large $B_0/V_0$, the tripole satellites lose their coherence and are not able to maintain large enough a shear on the vortex core and the whole structure tends towards axisymmetry.

4.2.2. Influence of stratification and of bottom topography

To assess the influence of stratification, we run simulations of the case $B0/V0=0.1,\alpha =2,l=2$, increasing $1/R_d^2$ from zero. We have run the model with three values of $1/R_d^2$: 0.25, 1.0, 4.0.

Figure 8 shows the various final states of vorticity and buoyancy (with increasing $1/R_d^2$ from left to right). Clearly, the vorticity satellites (in the periphery), contributing to the tripolar shape of the vortex, and the small-scale features are progressively erased by the increasing vortex stretching. The vortex core is also more circular. This confirms the linear stability analysis which indicated a decreasing instability with increasing $1/R_d^2$.

To assess the role of bottom topography, we run a simulation with a cubic exponential vortex ($\alpha =3.0$, with $B0/V0=0.1$, $1/R_d^2=0$, an elliptical perturbation ($l=2$) and a Gaussian topography $h=h_{0}exp(-{(r/r_0)}^2$, with $h_0=1.0,r_0=1.0$. In Figure 8, we compare the nonlinear evolution of a quartic exponential vortex ($\alpha =4.0$) perturbed elliptically ($l=2$) with $B_0/V_0=0.3$, over flat bottom (no topography) and over a seamount (positive topography). Over the flat bottom, the unstable vortex transforms first into a tripole with filamentary satellite, an evolution already seen (above). Finally, the vortex periphery becomes turbulent with many small-scale features. On the contrary, over the seamount the perturbation seems damped and the vortex becomes quasi circular. This stabilization is in agreement with the linear stability analysis.

Figure 9. End-states of simulations for a quartic exponential vortex ($\alpha =4.0$) with $B0/V0=0.3$ and $l=2$; (a) Vorticity and buoyancy maps at $t=0,100,200$ without bottom topography; (b) vorticity and buoyancy maps at $t=0,100,200$ over a Gaussian seamount of height $h_0=1.0$ and of radius $r_0=1.0$).

Figure 9. End-states of simulations for a quartic exponential vortex ($\alpha =4.0$) with $B0/V0=0.3$ and $l=2$; (a) Vorticity and buoyancy maps at $t=0,100,200$ without bottom topography; (b) vorticity and buoyancy maps at $t=0,100,200$ over a Gaussian seamount of height $h_0=1.0$ and of radius $r_0=1.0$).

4.2.3. Nonlinear evolutions with higher wavenumber perturbations

As shown by [28,29,30], more complex vortex aggregates than tripoles can exist in two-dimensional incompressible flows. Therefore, we study here their possible counterparts in the TQG model. We search quadrupolar vortices formed from the instability of circular vortices, with a mode $l=3$ initial perturbation. The growth of higher wavenumber perturbations on a power-exponential vortex, requires a higher degree of the power (a larger horizontal velocity shear), than that of low-mode perturbations. Here, we consider power-exponential vortices (equation 3) with $\alpha \ge 3.5$.

Figure 10a sumarizes the various nonlinear evolutions of the unstable vortices when $B_0/V_0=0.1$ or $0.3$. We note that a smaller horizontal shear is needed to stabilize a quadrupole nonlinearly with $B_0/V_0$ is small while quadrupoles degenerate into tripoles for larger values of $\alpha$ when $B_0/V_0=0.3$.

Again this shows that the presence of buoyancy in the mean flow weakens the ability of the initial perturbation to grow and (possibly) stabilize nonlinearly at large amplitude. For this larger value of $B_0/V_0$ the mode $l=3$ is less prevalent. It either decays and the vortex axisymmetrizes or it gives way to an elliptical (asymmetric) disturbance.

We also note that for large $B_0/V_0$ more filaments (short-scale perturbations) develop.

5. Discussion

Our study has, first, provided stability criteria and simple solutions to understand the effect of buoyancy coupling with vorticity, on the instability of circular vortices. They have shown that the vortices are then more unstable than in a 2D model, with notably larger growth rates for short waves. Second, and paradoxically, the numerical simulations have shown that, if the short waves indeed grew more in a TQG model, they have a rather stabilizing effect on the nonlinear evolution of the linearly unstable vortices, compared with the 2D case. This is explained as follows:

During the development of the instability of these vortices, the external annulus of vorticity breaks into l independent vortices. These peripheral vortices can deform and break the vortex core, leading either to rotating multipoles, or to l dipoles. But, in the TQG model, the appearance of filaments and of small vortices render these peripheral vortices less coherent. As a consequence, these latter exert a weaker shear on the vortex core and do not distort it as much. Thus the vortex is finally less unstable in the nonlinear evolution.

Note that these results do not depend on our explicit choice of alpha-exponential vortices for the mean circular flow; this choice is usual in view of observations of vortices at sea [9]. But previous studies [24] have shown that the analytical form of this radial profile of vorticity is not essential, as long as it belongs to a family of vortices (e.g. unshielded monopoles, monopoles, unshielded annuli, etc.).

Now, only barotropic (horizontal shear) instability has been studied here and it would be of interest to generalize our results to at least a two-layer TQG model to allow the study of baroclinic instability. Note that the TQG equations are not fully coupled (vorticity is not included in the buoyancy equation), contrary to the two-layer quasi-geostrophic (QG) equations in which each layerwise potential vorticity depends on the other layer streamfunction. Therefore, we can expect new coupling mechanisms between vorticity and buoyancy in a two-layer model.

Also, concerning the equations, the vorticity and buoyancy equations are different in mathematical structure (even if they were not coupled): buoyancy is directly advected while streamfunction is modified via vorticity advection followed by inversion; the inversion of the Laplacian has a widening effect on horizontal fields. This is another difference between the two-equation TQG model and the two-equation two-layer QG model.

Continuing with the analogy and differences between TQG and QG models, we must note that TQG dynamics share with SQG (surface quasi-geostrophy, [31]) the ability to generate small-scale features in the nonlinear dynamics. But a major dynamical difference exists in these two models. In the TQG model, the streamfunction derives from the inversion of the vorticity. Therefore, the external velocity field of a vortex decays as $1/r$. In the SQG model, the streamfunction is obtained from the inversion of the buoyancy (a square root of the Laplacian of streamfunction). This leads to a different Green’s function for the inversion and the external velocity field of a buoyancy patch (in SQG) decays in $1/r^2$. Therefore, the small-scale features are dynamically less active in the SQG than in the TQG model. Consequently, vortex instability in the SQG model has much similarity with that in the 2D model, contrary to TQG.

One can add that even in the shallow-water (SW) model, where velocity is not geostrophic, vortex instability is comparable to that observed in QG. Therefore, TQG stands apart in this whole range of models (2D, QG, SQG, SW, TQG). It must be remembered that it is an idealized and simplified model. Comparison with oceanic features is difficult, all the more so as the specificity of TQG is the formation of small-scale features, not easily observed via in situ nor satellite sensors.

6. Conclusions

The previous study of vortices in the thermal rotating shallow-water (TRSW) model [21] has yielded features and evolutions quite comparable to those observed in our study. A detailed comparison between TQG and TRSW is now desirable to assess the role of non-geostrophic velocity in vortex evolution in these models. Can our stability criteria be extended, for instance, to weakly non-geostrophic flows ?

For oceanographic applications, it could be of interest to compare vortex dynamics in the TQG model to those in a primitive-equation model with a surface mixed layer. Indeed, such a model would include all physical processes at work in the ocean: in particular low and high frequency dynamics (the latter being excluded in TQG), mechanical and thermohaline forcing by the atmosphere, complete thermodynamics, etc. Realistic ocean models could serve more easily a testbed for TQG than oceanographic observations (for the reasons mentioned above). A next step in the assessment of TQG models could be the coupling of two such models, one for the ocean surface and one for the atmospheric boundary layer, and the study of vortex and jet dynamics in such coupled models. This research is underway.