# Eady baroclinic instability of a circular vortex

Armand Vic, Xavier Carton, Jonathan Gula

► To cite this version:

Armand Vic, Xavier Carton, Jonathan Gula. Eady baroclinic instability of a circular vortex. Symmetry, 2022, 14 (7), pp.1438. 10.3390/sym14071438. hal-03926089

HAL Id: hal-03926089

https://cnrs.hal.science/hal-03926089

Submitted on 6 Jan 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

EADY BAROCLINIC INSTABILITY OF A CIRCULAR VORTEX

Armand Vic<sup>1,*</sup> ID<0000-0003-2046-5297, Xavier Carton<sup>1,*</sup> ID<0000-0002-7849-6611 and Jonathan Gula<sup>1,2,*</sup> ID<0000-0002-0876-9557

**Abstract:** The stability of two superposed buoyancy vortices is studied linearly in a two-level SQG model. The basic state is chosen as two top-hat vortices (with uniform buoyancy), coaxial and with same radius. Only the vertical distance between the two levels and the top and bottom buoyancy intensities are varied, the other parameters are fixed. The linear perturbation equations around this basic state form a two-dimensional ODE for which the normal and singular mode solutions are numerically computed. For normal modes, the system is stable if the vortices are sufficiently far from the other to prevent vertical interactions of the buoyancy patches, or if they are close to each other but with very different intensities, again preventing the resonance of Rossby waves around their contours. The vortex is unstable if the intensities are similar and if the vortices are close to each other vertically. The growth rates of the normal modes increase with the angular wave-number, also corresponding to shorter vertical distances. The growth rates of the singular modes do not depend much on the bottom buoyancy at short time, but, as expected, they converge towards the growth rates of the normal modes. This study remaining linear does not predict the final evolution of such unstable vortices. This nonlinear evolution will be studied in a sequel of this work.

**Keywords:** SQG; Rankine vortices; normal modes; singular modes; linearisation around basic state.

# 1. Introduction

Vortices are energetic features in many turbulent flows. Paramount among them are geophysical flows, where vortices play an essential role in the planetary transport of energy, heat, moisture for atmospheric vortices, and salinity for oceanic vortices. Vortices are long-lived recirculation motion with a lifetime longer than their turnover period. Therefore it is essential to study the mechanisms underlying their robustness, or their possible destabilization. This has been the subject of many papers in the past ([1–5] and references therein). But many of these papers considered vortex instability in a layered model of the ocean (often corresponding to the ocean above the main thermocline, i.e. above 500 m depth, and below it). This problem is then called the Phillips baroclinic instability of these vortices. The focus was then on fairly large vortices in the ocean (vortices wider than 30 km in radius). Far fewer papers were devoted to the study of vortex stability in a level model of the ocean, where only the density interfaces are concerned. This problem is called the Eady baroclinic instability of these vortices. Such interfaces (the ocean surface, the ocean bottom or the thermocline) play an essential role in ocean dynamics. It was shown recently that a model describing these surfaces only, can represent the dynamics of smaller vortices (with radii 10-30 km) which are abundant in the ocean. Such a model is the Surface Quasi-Geostrophic (SQG) model, employed in this study. This model describes the time evolution of buoyancy anomalies on surfaces, in a rotating stratified flow, with null internal potential vorticity. The 3D internal dynam-

**Citation:** Vic, A.; Carton, X.; Gula, J. Eady vortex problem. *Symmetry* **2021**, *1*, 0. https://doi.org/

Received: XXX

Accepted: XXX

Published: XXX

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Submitted to *Symmetry* for possible open access publication under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ UBO, CNRS, Ifremer, IRD, Laboratoire d'Océanographie Physique et Spatiale (LOPS), IUEM, 29280 Plouzané, France; lops-dir@listes.ifremer.fr</footnote>

<footnote n="2">² Institut Universitaire de France (IUF), Paris, France; secretariat.iuf@recherche.gouv.fr</footnote>

<footnote>* Correspondence: armand.vic@univ-brest.fr (A.V.) ; xcarton@univ-brest.fr (X.C. ) ; jonathan.gula@univ-brest.fr (J.G.)</footnote>
ics (vertically, between the horizontal surfaces) are driven by the buoyancy anomalies
on these surfaces.

The previous studies of vortex stability in the SQG model concerned the horizontal
shear (barotropic) instability of a single vortex in one- or two-level configurations. The
analytical stability of two superposed vortices in a two-level SQG model has not been
investigated before. Badin and Poulin [6] or Harvey and Ambaum [7] studied the
barotropic instability of a single vortex. Here, we study the baroclinic instability (ver-
tical shear instability of a rotating fluid) of two superimposed vortices. The two-level
SQG model (see [8]) is adapted to deal analytically with this vortex instability problem.

Analytically, the model solves a hyperbolic equation for the transport of buoyancy of the
surfaces and a 3D Laplacian equation on the streamfunction (an elliptical equation to
invert the buoyancy distribution into a flow field) with Neumann boundary condition.

The model equations and their numerical implementation are developed in section 2.

Section 3 presents the basic state composed of two vortices (one at the ocean surface, one at the ocean thermocline or bottom). These top-hat vortices can have different intensities and their vertical separation can be varied. In section 4, we linearize the equations around this steady state and we study how the perturbation grows. This perturbation can be a normal mode or a singular mode. Section 4 presents and interprets the results.

A conclusion and perspectives follow. Two appendices present details about analytical
and numerical computations.

## 2. Surface quasi-geostrophic model and equations

The framework is a two-level surface quasi-geostrophic (SQG) model. The two
horizontal surfaces are the actual surface and the bottom of the ocean (or the thermocline
if the temperature - and density - gradient is sufficiently abrupt at this depth [9,10]);
they are vertically separated by a height *H* (see Figure 1). The buoyancy (or potential
temperature) distributions are contained in the two levels. They are connected by a
condition of null potential vorticity inside the ocean.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Scheme of the two horizontal layers.</figcaption></figure>

The surface quasi-geostrophic (SQG) model is the restriction of the complete quasi-
geostrophic model – introduced by Charney in 1948[11] – to null internal potential
vorticity distributions (for more details on the SQG model, see [8,12]). Potential vorticity
is thus concentrated as a vertical Dirac distribution, reducing to a planar buoyancy
anomaly, at the two (upper and lower) boundaries.

Assuming constant Brunt-Vaïsala and Coriolis frequencies, the quasi-geostrophic model is governed by the conservation of potential vorticity $q$ in the fluid volume, in the absence of forcing and of dissipation for the flow:

$$dq/dt = \partial_t q + J(\psi, q) = 0, \quad \text{for } 0 < z < H,$$

associated with following 3D Laplace equation:

$$q = \nabla_h^2 \psi + \frac{f_0^2}{N_0^2} \partial_z^2 \psi, \quad \text{for } 0 < z < H, \qquad (1)$$

where $J$ is the horizontal Jacobian operator and $\psi$ is the streamfunction (remember that the horizontal velocity is $u = -\partial_y \psi, v = \partial_x \psi$). The boundary conditions of the 3D model at the top and bottom of the domain, are the horizontal advection of the buoyancies $b$ on these surfaces:

$$[\partial_t + J(\psi, \cdot)]b = 0, \quad z = 0, H \qquad (2)$$

with $b = f_0 \partial_z \psi$.

The surface quasi-geostrophic equations are therefore the restriction of this model to $q = 0$ in the fluid interior. This leads to a model defined only in terms of the surface and bottom buoyancies, related to streamfunction as above.

From these equations, one can define a horizontal length scale, called the internal radius of deformation, via $\sigma = \frac{N_0 H}{f_0}$. This scale represents the distance over which the buoyancy and Coriolis accelerations have similar intensities. Now we normalize in [0, 1] the vertical scale ($H$) to have the following set of equations:

$$\left\{ \begin{array}{ll} \nabla_h^2 \psi + \frac{1}{\sigma^2} \partial_z^2 \psi = 0 & \text{for } 0 < z < 1 \\ \left. \frac{\partial \psi}{\partial z} \right|_{z=0} = b^s, & \frac{Db^s}{Dt} = 0 \\ \left. \frac{\partial \psi}{\partial z} \right|_{z=1} = b^b, & \frac{Db^b}{Dt} = 0 \end{array} \right. \qquad (3)$$

where $\frac{D}{Dt}$ is the horizontal Lagrangian derivative and the superscripts s and b represent respectively “surface” and “bottom”.

The first equation of system (3) in horizontal Fourier space ($k, l, z$) gives

$$\frac{\partial^2 \hat{\psi}}{\partial z^2} = K^2 \sigma^2 \hat{\psi} \qquad (4)$$

where $K^2 = k^2 + l^2$. This equation with boundary conditions for buoyancies gives in Fourier space

$$\hat{\psi}(k, l, z) = \frac{1}{K\sigma \sinh(K\sigma)} \left( \hat{b}^b \cosh(K\sigma z) - \hat{b}^s \cosh(K\sigma(1-z)) \right) \qquad (5)$$

The mean flow in our SQG model are two top-hat vortices (i.e. two vortices with constant buoyancy in a disk of radius unity). Here we analyse the linear stability of this mean flow (or basic state). This problem is called the Eady baroclinic instability of this vortex. With this mean flow geometry, the polar coordinates are a natural choice.

### 3. Mean flow calculation

Here we calculate the flow field associated with these two top-hat vortices. Firstly we remind the form of cylindrical Fourier transforms.

#### 3.1. Preliminaries about Fourier decomposition in cylindrical coordinates

In cylindrical coordinates, consider a function $f(r, \phi, z)$ sufficiently regular.
• Then $f$ is $2\pi$-periodic in $\phi$ so can be decomposed in Fourier modes :

$$f(r, \phi, z) = \sum_{n \in \mathbb{N}} \tilde{f}(r, n, z) e^{in\phi} \qquad (6)$$

with

$$\tilde{f}(r, n, z) = \frac{1}{2\pi} \int_{0}^{2\pi} f(r, \phi, z) e^{-in\phi} d\phi \qquad (7)$$

• For every $n \in \mathbb{N}$ and $z \in \mathbf{R}^+$, $\tilde{f}(\cdot, n, z)$ are functions which can be written as inverse Hankel transforms :

$$\tilde{f}(r, n, z) = \int_{0}^{\infty} \hat{f}(\rho, n, z) J_{n}(\rho r) \rho d\rho \qquad (8)$$

where $J_n$ are the Bessel functions and with

$$\hat{f}(\rho, n, z) = \int_{0}^{\infty} \tilde{f}(r, n, z) J_{n}(\rho r) r dr \qquad (9)$$

In fine, the function can be decomposed as

$$f(r, \phi, z) = \sum_{n \in \mathbb{N}} \int_{0}^{\infty} \hat{f}(\rho, n, z) J_{n}(\rho r) \rho d\rho e^{in\phi} \qquad (10)$$

with

$$\hat{f}(\rho, n, z) = \frac{1}{2\pi} \int_{0}^{2\pi} \int_{0}^{\infty} f(r, \phi, z) J_{n}(\rho r) r dr e^{-in\phi} d\phi \qquad (11)$$

88 where $\hat{f}$ is the Fourier transform of $f$ in cylindrical Fourier coordinates.

89 **Remark 1.** For the second bullet point, the Bessel function could be arbitrary but the n-th function is retained because it is a solution of the Laplace equation.

90 3.2. Application to SQG flows

91 We decompose $\psi$ as in (10) to get :

$$\psi(r, z, \phi) = \sum_{n \in \mathbb{N}} \int_{0}^{\infty} \hat{\psi}(\rho, n, z) J_{n}(\rho r) \rho d\rho e^{in\phi} \qquad (12)$$

92 with $\hat{\psi}$ the horizontal Fourier transform of $\psi$, already computed in (5). We deduce :

$$\psi = \sum_{n \in \mathbb{N}} \int_{0}^{\infty} \frac{\rho J_{n}(\rho r)}{\rho \sigma \sinh(\rho \sigma)} (\hat{b}^{\mathrm{b}}(\rho, n, t) \cosh(\rho \sigma z) - \hat{b}^{\mathrm{s}}(\rho, n, t) \cosh(\rho \sigma(1-z))) d\rho e^{in\phi}. \qquad (13)$$

93 So at the two boundaries (surface and bottom) we have :

$$\psi^{\mathrm{s}}(r, \phi, z = 0, t) = \sum_{n \in \mathbb{N}} \int_{0}^{\infty} \frac{J_{n}(\rho r)}{\sigma \sinh(\rho \sigma)} (\hat{b}^{\mathrm{b}} - \hat{b}^{\mathrm{s}} \cosh(\rho \sigma)) d\rho e^{in\phi} \qquad (14a)$$

$$\psi^{\mathrm{b}}(r, \phi, z = 1, t) = \sum_{n \in \mathbb{N}} \int_{0}^{\infty} \frac{J_{n}(\rho r)}{\sigma \sinh(\rho \sigma)} (\hat{b}^{\mathrm{b}} \cosh(\rho \sigma) - \hat{b}^{\mathrm{s}}) d\rho e^{in\phi} \qquad (14b)$$

94 **Remark 2.** From now, we will denote by capital letters the basic state variables and by lowercase letters the perturbed variables. For example, the total streamfunction at the surface will be

95 $\Psi^s + \psi^s$.
3.3. Basic state : two top-hat vortices

We take as basic states two top-hat vortices (i.e two circular plateaus of constant buoyancy) with the same dimensionless radius $R = 1$ but not the same intensity.

$$B^s = B_0^s(1 - H(r-1)) \qquad (15a)$$

$$B^b = B_0^b(1 - H(r-1)) \qquad (15b)$$

where $H(x) = \begin{cases} 0 & \text{if } x < 0 \\ 1 & \text{if } x > 0 \end{cases}$ is the Heavyside function.

To have the streamfunction of the basic state, we need the Fourier transforms of the buoyancies : for $i = s, b$, we have $\hat{B}^i(\rho, n) = 0$ for $n \neq 0$ because $B^i$ is independent of $\phi$. For $n = 0$, because $\frac{d}{dx}(xJ_1(x)) = xJ_0(x)$, we have $\hat{B}^i(\rho, 0) = B_0^i \frac{J_1(\rho)}{\rho}$. Then the streamfunction for the basic state at the two boundary levels is :

$$\Psi^s = \int_0^\infty \frac{J_1(\rho)J_0(\rho r)}{\rho\sigma} \left( \frac{B_0^b}{\sinh(\rho\sigma)} - \frac{B_0^s}{\tanh(\rho\sigma)} \right) d\rho \qquad (16a)$$

$$\Psi^b = \int_0^\infty \frac{J_1(\rho)J_0(\rho r)}{\rho\sigma} \left( \frac{B_0^b}{\tanh(\rho\sigma)} - \frac{B_0^s}{\sinh(\rho\sigma)} \right) d\rho \qquad (16b)$$

The flow at the surface of the ocean induced by the vortices is along $\vec{e}_\phi$ and because $J_0' = -J_1$, we have :

$$U_r^s = -\frac{1}{r} \partial_\phi \Psi^s = 0 \qquad (17a)$$

$$U_\phi^s = \partial_r \Psi^s = \int_0^\infty \frac{J_1(\rho)J_1(\rho r)}{\sigma} \left( \frac{B_0^s}{\tanh(\rho\sigma)} - \frac{B_0^b}{\sinh(\rho\sigma)} \right) d\rho \qquad (17b)$$

Similarly, the flow at the bottom of the ocean is :

$$U_r^b = 0 \qquad (18a)$$

$$U_\phi^b = \int_0^\infty \frac{J_1(\rho)J_1(\rho r)}{\sigma} \left( \frac{B_0^s}{\sinh(\rho\sigma)} - \frac{B_0^b}{\tanh(\rho\sigma)} \right) d\rho \qquad (18b)$$

We now introduce the quantities :

$$I_n(r, \sigma) := \int_0^\infty \frac{J_n(\rho)J_n(\rho r)}{\sigma \tanh(\rho\sigma)} d\rho \qquad (19a)$$

$$M_n(r, \sigma) := \int_0^\infty \frac{J_n(\rho)J_n(\rho r)}{\sigma \sinh(\rho\sigma)} d\rho \qquad (19b)$$

such that

$$U_\phi^s = B_0^s I_1 - B_0^b M_1 \qquad (20a)$$

$$U_\phi^b = B_0^s M_1 - B_0^b I_1 \qquad (20b)$$

**Remark 3.** This is indeed a steady basic state : since there is no radial velocity, the buoyancy anomaly, which is a tracer, will remain a circular patch if unperturbed.
### 4. Evolution of the vortex boundaries in the linear instability of the vortex

Now we perturb each vortex by deforming its contour. Because $\frac{D_b^{s,b}}{Dt} = 0$, an initial plateau in buoyancy will remain such at all times (with a deformed external contour).

Therefore the lateral jump in buoyancy will always exist and we can define the vortex boundaries as the place where the jump lies. The evolution of the vortex boundaries will measure the stability of this particular basic state.

Assume that the radii $R^s$ and $R^b$ of the vortices at the surface and the bottom are disturbed from their basic states 1, as represented in Figure 2.

**Remark 4.** We assume that during the linear stage of instability, as the perturbation amplitude remains small, the boundary will not be locally multi-valued, such that we can use the following parameterisation:

where $\eta^i$ is small compared with 1.

$$ \begin{cases} R^s(\phi, t) = 1 + \eta^s(\phi, t) \\ R^b(\phi, t) = 1 + \eta^b(\phi, t) \end{cases} \qquad (21) $$

<figure><img src="image_3.png" /><figcaption>Figure 2. Perturbation of the buoyancy disk, with $\eta \ll 1$.</figcaption></figure>

Then the buoyancy at the surface is

$$ B_0^s(1 - H(r - R^s(\phi, t))) = \begin{cases} B_0^s & \text{if } r < R^s \\ 0 & \text{if } r > R^s \end{cases} \qquad (22) $$

where $H$ is the Heavyside function. A similar form can be derived for the buoyancy

at the bottom.

Because we chose (15) as basic buoyancies, the perturbed buoyancy at the surface is

130:

$$ b^s = B_0^s(H(r-1) - H(r-1-\eta^s)) = \begin{cases} B_0^s & \text{if } 1 < r < 1+\eta^s \\ -B_0^s & \text{if } 1+\eta^s < r < 1 \\ 0 & \text{else} \end{cases} \qquad (23) $$

As a distribution, for small $\eta^s$, we have
$$ \begin{cases} b^s(r, \phi, t) = B_0^s \eta^s(\phi, t) \delta_1(r) \\ b^b(r, \phi, t) = B_0^b \eta^b(\phi, t) \delta_1(r) \end{cases} \qquad (24) $$

132 where $\delta_1$ is the Dirac distribution in 1.

133 Because the perturbed streamfunctions from Eq. (13) are searched, we need the 134 Fourier transform of the perturbed buoyancies:

$$ \hat{b}^s(\rho, n, t) = \frac{1}{2\pi} \int_0^{2\pi} \int_0^{\infty} B_0^s \eta^s(\phi, t) \delta_1(r) J_n(\rho r) r dr e^{-in\phi} d\phi \quad (25) $$

$$ = \frac{B_0^s}{2\pi} \int_0^{2\pi} J_n(\rho) \eta^s(\phi, t) e^{-in\phi} d\phi \qquad (26) $$

$$ \hat{b}^s(\rho, n, t) = B_0^s J_n(\rho) \hat{\eta}^s(n, t) \qquad (27) $$

135 and a similar formula for $\hat{b}^b(\rho, n, t)$. Because the equations are linear, the perturbed 136 streamfunctions can be computed from the formula (14):

$$ \psi^s(r, \phi, t) = \sum_{n \in \mathbb{N}} \int_0^\infty \frac{J_n(\rho r)}{\sigma \sinh(\rho \sigma)} (\hat{b}^b - \hat{b}^s \cosh(\rho \sigma)) d\rho e^{in\phi} \quad (28) $$

$$ \psi^s(r, \phi, t) = \sum_{n \in \mathbb{N}} (B_0^b \hat{\eta}^b M_n - B_0^s \hat{\eta}^s I_n) e^{in\phi} \quad (29) $$

$$ \psi^b(r, \phi, t) = \sum_{n \in \mathbb{N}} (B_0^b \hat{\eta}^b I_n - B_0^s \hat{\eta}^s M_n) e^{in\phi} \quad (30) $$

137 where $I_n(r, \sigma)$ and $M_n(r, \sigma)$ are defined in (19).

138 The (total and perturbed) radial flow at the boundary of the surface vortex is on the 139 one hand:

$$ u_r^s(R^s(\phi, t), \phi, t) = - \frac{1}{1 + \eta^s} \partial_\phi \psi^s (1 + \eta^s, \phi, t) \quad (31) $$

141 and on the other hand

$$ u_r^s(R^s(\phi, t), \phi, t) = \frac{DR^s}{Dt} = \frac{\partial R^s}{\partial t} + \frac{1}{R^s} (U_\phi^s + u_\phi^s) \frac{\partial R^s}{\partial \phi} \quad (32) $$

142 The justification of the first equality $u_r^s(R^s(\phi, t), \phi, t) = \frac{DR^s}{Dt}$ is the following: $R^s$ is the radial coordinates of a material line (the boundary of the vortex). So its rate of Lagrangian displacement in time is exactly the radial velocity of the flow.

143 With the equality of Equations (31) and (32), using a Taylor expansion for small 144 amplitude perturbations (details are given in Appendix A), we obtain the following 145 system:

$$ \begin{cases} \partial_t \eta^s = -\partial_\phi \psi^s - U_\phi^s \partial_\phi \eta^s \\ \partial_t \eta^b = -\partial_\phi \psi^b - U_\phi^b \partial_\phi \eta^b \end{cases} \qquad (33) $$

149 where the functions are applied in $r = 1$.

150 Assuming that one mode of perturbation will grow faster than any other, we retain 151 only one Fourier mode $\eta^{s,b}(\phi, t) = \hat{\eta}^{s,b}(t) e^{in\phi}$ so that we obtain the matrix form with 152 $\hat{\eta} = (\hat{\eta}^s \quad \hat{\eta}^b)$:

153
$$ \partial_t \hat{\eta} = in \begin{pmatrix} B_0^s(I_n - I_1) + B_0^b M_1 & -B_0^b M_n \\ B_0^s M_n & B_0^b(I_1 - I_n) - B_0^s M_1 \end{pmatrix} \hat{\eta} \quad (34) $$

## 5. Results

### 5.1. Preliminaries: study of the integrals $I_n - I_1$ and $M_n$

We failed to compute analytically the integrals $I_n - I_1$ and $M_n$ so from now, we do a numerical study of the stability. Numerically, these integrals have poor convergence properties. In particular the integral $I_n - I_1$ is not absolutely convergent. We present the method to compute it in Appendix B. Nevertheless, we obtain a graphical representation of $I_n - I_1$ and $M_n$ with respect to $\sigma$ in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3. Graphs of the integrals $I_2 - I_1$ and $M_2$ with respect to $\sigma$.</figcaption></figure>

### 5.2. Normal modes

In this section, we consider normal mode perturbations to the vortex boundary. This means that the time dependence of the $\eta^{s,b}$ is $\hat{\eta}^{s,b}(t) = \mu^{s,b} e^{-i\omega_n t}$ with $\omega_n = a_n + ib_n \in \mathbb{C}$ where $b_n$ is the growth rate and $\mu^{s,b} \in \mathbb{R}^+$. In order to conclude about stability, we are now interested in the sign of $b_n$, the imaginary part of $\omega_n$. Thanks to (34), we obtain an eigenvalue problem $A_n\mu = -\frac{\omega_n}{n}\mu$ where $\mu = \begin{pmatrix} \mu^s \\ \mu^b \end{pmatrix}$ and

$$ A_n = \begin{pmatrix} B_0^s(I_n - I_1) + B_0^b M_1 & -B_0^b M_n \\ B_0^s M_n & B_0^b(I_1 - I_n) - B_0^s M_1 \end{pmatrix}. \quad (35) $$

Because $A_n \in \mathcal{M}_n(\mathbb{R})$, there are two possibilities for the eigenvalues. They can be real, then $b_n = 0$ and the basic state has a neutral stability; or they can be complex conjugate, then one of the two eigenvalues has $b_n > 0$ and the basic state is unstable.
**Remark 5.** The normal mode $n = 1$ is always stable because

$$A_1 = \begin{pmatrix} B_0^b M_1 & -B_0^b M_1 \\ B_0^s M_1 & -B_0^s M_1 \end{pmatrix} \qquad (36)$$

Conclusions on this flow stability are obtained by computing the discriminant of $\chi_n(X)$, the characteristic polynomial of $A_n$:

$$\chi_n(X) = X^2 + (B_0^b - B_0^s)(I_n - I_1 - M_1)X \qquad (37)$$

$$-B_0^b B_0^s ((I_n - I_1)^2 + M_1^2 - M_n^2) - (I_n - I_1)M_1 ((B_0^s)^2 + (B_0^b)^2). \qquad (38)$$

To know if the roots of this polynomial are real or complex conjugate, we compute the sign of the discriminant $\Delta$:

$$\Delta = (B_0^b - B_0^s)^2 (I_n - I_1 - M_1)^2 \qquad (39)$$

$$+ 4B_0^b B_0^s ((I_n - I_1)^2 + M_1^2 - M_n^2) \qquad (40)$$

$$+ 4(I_n - I_1)M_1 ((B_0^s)^2 + (B_0^b)^2) \qquad (41)$$

$$\Delta = (I_n - I_1 + M_1)^2 (B_0^b + B_0^s)^2 - 4B_0^b B_0^s M_n^2 \qquad (42)$$

<figure><img src="image_5.png" /><figcaption>Figure 4. $b_n$ for $n = 2, 3, 4$ and 5, with respect to $\sigma$ and $\frac{B_0^s}{B_0^b}$.</figcaption></figure>

<footnote n="166">166 has two real eigenvalues : 0 and $M_1(B_0^b - B_0^s)$.</footnote>

<footnote n="167">167 The conclusion is : if $\Delta > 0$, then the system is neutral because $b_n = 0$ and if $\Delta < 0$, then the growth rate $b_n = \frac{n\sqrt{-\Delta}}{2} \neq 0$ and the system is unstable.</footnote>

<footnote n="169">169 **Remark 6.** If we take, as Badin and Poulin in [6] $B_0^b = 0$, we obtain $\Delta = B_0^s^2(I_n - I_1)^2 > 0$ and then we recover their dispersion relation $\omega_n = nB_0^s(I_1 - I_n) \in \mathbf{R}$.</footnote>

<footnote n="171">171 Because the case $n = 1$ is always stable, we plot in Figure 4 the normal modes for $n = 2, 3, 4$ and 5.</footnote>
The dark purple zone is where $b_n = 0$. There, $\Delta > 0$ and the system reaches a stable state with the following dispersion relation:

$$ \omega_n(\sigma) = \frac{n}{2} \left[ (B_0^s - B_0^b)(I_n - I_1 - M_1) \pm \sqrt{\Delta} \right] \quad (43) $$

The four angular modes have similar stability properties on the top-hat vortex, but for different values of the physical parameters ($\sigma, B_0^s/B_0^b$). For each mode, we can separate three stable zones:

* in the right hand side of each panel, where $\sigma$ is larger than a threshold $\sigma_{critic}$ depending on $n$. We recover here the results of [6] or [7]. They found that a top-hat vortex, alone in a SQG model, is stable. In this area, the system is linearly stable for barotropic (horizontal shear) instability. They are sufficiently far from the other ($\sigma$ is proportional to $H$) so we could neglect the interactions. The two-layer SQG model is then viewed as two one-layer SQG models where there are two independent top-hat vortices. We define $\sigma_{critic}$ as the critical value of $\sigma$ leading to instability, all other parameters being fixed. $\sigma_{critic}$ is a decreasing function with respect to $n$. The stability of high mode perturbations is reached for a smaller distance between vortices than low mode perturbations. This is due to the relation between horizontal and vertical wave numbers in the SQG model.

* in the top left and the bottom left sides of each panel, the system is also stable. In these areas, the vortices are close to each other but have very different intensities. An interpretation of this situation could be that perturbations on one of the vortice have very different phase speeds around the contour than for the other vortex. The impossibility for these two (Rossby) waves to phase lock and resonate stabilizes the whole system.

The system is unstable if the mean buoyancy intensities are similar and if the vortices are vertically close to each other.

For a given fixed mode and a fixed ratio $B_0^s/B_0^b > 1$, an interpretation of the change of stability is the following:

* for small $\sigma$, the two vortices are too close to each other for the wave to grow; thus short wave cut off (usual for the Eady model) can be explained by the absence of phase locking between waves.

* for intermediate $\sigma$, the distance between the two vortices allows the mode to grow (phase locking with the proper phase shift is possible) and then the system is unstable. The smaller $\sigma$ is, the shorter the most unstable waves are.

* for large $\sigma$, the two vortices are far from each other, wave-wave interaction is weak and the mode is stable.

### 5.3. Singular modes

System (34) is a 2 × 2 system, and the matrix $A_n$ is independent of time $t$. The solution is then given by $\hat{\sigma} = e^{A_n t} \hat{\sigma}_0$. Since matrix $A$ is not self-adjoint, linear combinations of normal modes can grow faster than the normal modes [13]; they are the singular modes of the problem. Therefore we calculate the singular modes of stability of this vortex flow. They are defined as the maximal growth rate of a given norm of these perturbations. Here we choose the squared perturbation buoyancy as the norm. The numerical method for singular mode calculation given in the appendix of [14] is implemented here. The growth rates of the singular modes are the eigenvalues of the matrix $e^{A_n t}$. These growth rates are shown with respect to the same physical parameters ($\sigma, B_0^s/B_0^b$) as the normal modes, for increasing values of the time $t$ and for $n = 2$ on figure 5.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Singular mode for $n = 2$ for different times $t = 1, 2, 5, 10, 50$. Here $B_0^s = 0.25$. The bottom right panel represents the growth rates for normal mode $n = 2$. We can note the convergence of singular modes to normal modes.</figcaption></figure>

For short time, the singular modes growth rates are concentrated near the region of small $\sigma$. Indeed, at small times, only short waves can grow corresponding to buoyancy surfaces close to each other vertically. Also, for small $\sigma$, the singular mode instability occurs for every $B_0^b$, i.e. even a weak bottom buoyancy is sufficient then to allow the phase locking of these short waves. Also this implies that two vertically close vortices, with different mean buoyancies, are unstable for singular modes but stable in normal modes for small time. Note that a similar remark was made in [15] about the independence of singular mode growth rates to the barotropic component of the flow at short times. As time grows, the singular mode growth rates for $n = 2$ converge towards those of the corresponding normal modes; this result is valid for the other $n \in \mathbb{N}$. This confirms the result of previous studies [13–15].
## 6. Conclusions and perspectives

In this study, we have developed analytically and numerically the calculation of growth rates for the instability of two superposed vortices. The theory and the computation of the mean flow were first done analytically. The computation of the normal and singular modes was done numerically. The two-level SQG model and the considered steady state are idealized, but the provide simple results for the Eady baroclinic instability of two superimposed vortices: stability for vertically distant vortices, instability for vertically close vortices, similar in intensity and instability in singular modes only for small time for close and different in intensity vortices.

Though these results pertain to an idealized vortex, we can apply them to the ocean. Using the following values $f_0 = 10^{-4}s^{-1}$, $N_0 = 5 \ 10^{-3}s^{-1}$, $R = 2.5 \ 10^4m$, $h = 10^3m$, $V = 0.5m/s$, where V is the rotational velocity of an oceanic vortex, we obtain the following length and time scales, $L = 2.5 \ 10^4m$, $T = 6.3 \ 10^4s$. Firstly, we can use the values of $\sigma$ to determine which vortices can be unstable: strong instability occurs for $\sigma = 0.2$ leading to most unstable vortices having a thickness of 100 m which indeed corresponds to small vortices (with radii close to 10 km). Secondly, we can compute the growth rates of such normal mode perturbations in the ocean: in dimensionless terms, they are on the order of 0.3 for $n = 2$. This corresponds to a typical time scale for the growth of these perturbations of $6.3/0.3 \ 10^4s$, about 2.5 days. This timescale is slightly shorter than that found in the two-layer Phillips problem of mesoscale vortex baroclinic instability, about 4 days [1,15].

Finally, we must also note that we have studied only the linear instability of such vortices. The natural follow-up of these calculations is the study of the long-term, nonlinear evolution of these unstable vortices. This will indicate if the linearly unstable waves found here, can be stabilised on the long run via nonlinear wave interactions and which shape the nonlinearly stabilised vortices would take. Furthermore, this will allow the variation of several parameters not included here:

* investigate the effect of different radii for the two vortices

* shift one vortex with respect to the other and study the evolution of tilted vortices

* consider two different modes $n_s$ and $n_b$ of perturbation for the two vortices

* consider other radial shapes for the vortices (Gaussian,...)

This numerical study is under way and will complement this first paper.

**Funding:** This research has been supported by external funding: the first author has received a Ph.D grant from Ecole Normale Superieure, Rennes

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

## Abbreviations

The following abbreviation is used in this manuscript:

SQG Surface Quasi Geostrophic

## Appendix A. Proof of the system (33)

We will develop the calculus for the surface only. The bottom case is similar. On one hand, from Equation (31), we have:

$$u_r^s(R^s(\phi, t), \phi, t) = -(1 + O(\epsilon)) (\partial_\phi \psi^s(1, \phi, t) + O(\epsilon^2)) \quad (A1)$$

$$= -\partial_\phi \psi^s(1, \phi, t) + O(\epsilon^2) \quad (A2)$$
On the other hand, from Equation (32) we have :

$$u_r^s(R^s(\phi, t), \phi, t) = \frac{\partial \eta^s}{\partial t}(\phi, t) + \frac{1}{1 + \eta^s(\phi, t)} \quad (A3)$$

$$[U_{\phi}^{s}(1 + \eta^{s}(\phi, t), \phi, t) + U_{\phi}^{s}(1 + \eta^{s}(\phi, t), \phi, t)] \frac{\partial \eta^{s}}{\partial \phi}(\phi, t) \quad (A4)$$

$$= \partial_t \eta^s(\phi, t) + (1 + O(\epsilon)) [U_{\phi}^s(1, \phi, t) + O(\epsilon)] \partial_{\phi} \eta^s(\phi, t) \quad (A5)$$

$$= \partial_t \eta^s(\phi, t) + U_{\phi}^s(1, \phi, t) \partial_{\phi} \eta^s(\phi, t) + O(\epsilon^2) \quad (A6)$$

And then, at the order $O(\epsilon)$, we obtain

$$\partial_t \eta^s(\phi, t) = -\partial_\phi \psi^s(1, \phi, t) - U_\phi^s(1, \phi, t) \partial_\phi \eta^s(\phi, t) \quad (A7)$$

**Remark A1.** Note a difficulty we did not mention earlier : $U_{\phi}^s$ is not differentiable in the classical way. Formally, we should work with smooth approximation of top-hat vortices and then move on to the limit.

**Appendix B. Proof of convergence and numerical method to compute $I_n - I_1$ and $M_n$**

Recall $I_n - I_1 = \frac{1}{\sigma} \int_0^\infty f_n(x)dx$ and $M_n = \frac{1}{\sigma} \int_0^\infty g_n(x)dx$ with $f_n(x) = \frac{I_n(x)^2 - I_1(x)^2}{\tanh(\sigma x)}$ and $g_n(x) = \frac{I_n(x)^2}{\sinh(\sigma x)}$ (see Figure A1). The two functions are continuous on $]0, \infty[$. Let's do the analysis of convergence in 0 and in $\infty$:

* For $x$ in a neighborhood of 0, $J_n(x) \sim \frac{x^n}{2^n n!}$ so for $n > 1$:

$$f_n(x) \sim \frac{J_1(x)^2}{\tanh(\sigma x)} \quad (A8)$$

$$f_n(x) \sim \frac{x}{4\sigma} \quad (A9)$$

and for $n \ge 1$:

$$g_n(x) \sim \frac{x^{2n-1}}{4^n (n!)^2 \sigma}. \quad (A10)$$

So the two functions are integrable in 0.

* For $x$ in a neighborhood of $+\infty$,

$$J_n(x) = \sqrt{\frac{2}{\pi x}} \sin\left(x - \frac{n\pi}{2} + \frac{\pi}{4}\right) - \frac{4n^2 - 1}{4\sqrt{2\pi x}^{\frac{3}{2}}} \sin\left(x - \frac{n\pi}{2} - \frac{\pi}{4}\right) + o\left(\frac{1}{x^{\frac{3}{2}}}\right) \quad (A11)$$

The computation gives for $n = 2p > 1$:

$$f_{2p}(x) \sim \frac{2 \sin(2x)}{\pi x} \quad (A12)$$

and for $n = 2p + 1 \ge 1$:

$$f_{2p+1}(x) \sim \frac{1 - (2p + 1)^2}{\pi x^2} \cos(2x). \quad (A13)$$

We can quickly conclude for the odd case because $f_{2p+1} = O(\frac{1}{x^2})$ is absolutely convergent in $+\infty$. The even case is a modified integral sine so that it converges. For every $n \in \mathbb{N}$, we have a quick convergence in $+\infty$ for $g_n$:

$$g_n(x) \sim \frac{4}{\pi x} \sin^2\left(x - \frac{n\pi}{2} + \frac{\pi}{4}\right) e^{-\sigma x}. \quad (A14)$$
<figure><img src="image_7.png" /><figcaption>**Figure A1.** The integrands $f_n$ and $g_n$ are in solid line, the asymptotes are plotted with crosses and the envelopes for the top right panel are in dashed lines. Notice that for the bottom right panel, the asymptote depends only on the parity of $n$. This explains why there are only two asymptotes plotted. Here we take the parameter $\sigma = 1$.</figcaption></figure>

Numerically, the only difficult point is the fact that the integral sine does not absolutely converge so that the classically implemented methods to compute integrals are not adapted. A python routine exists to compute integral sine and this is what we use. The idea is to cut the integrals in three parts, $\int_0^\infty = \int_0^\epsilon + \int_\epsilon^A + \int_A^\infty$, to use approximation for the integrals in 0 and in $+\infty$ and to use classical python routine in $[\epsilon, A]$. With the asymptotic developments we used, we obtain:

* for $f_n$ in 0:

$$ \int_0^\epsilon f_n(x)dx \approx -\frac{\epsilon^2}{8\sigma} \qquad (A15) $$

* for $f_{2p}$ in $+\infty$:

$$ \int_A^\infty f_{2p}(x)dx \approx \frac{2}{\pi} \int_{2A}^\infty \frac{\sin t}{t}dt \qquad (A16) $$

$$ \approx 1 - \frac{2}{\pi}\mathrm{Si}(2A) \qquad (A17) $$

<footnote n="284">284 where $\mathrm{Si}(x) = \int_0^x \frac{\sin(t)}{t}dt$ is the integral sine function.</footnote>
• for $f_{2p+1}$ in $+\infty$:

$$ \int_A^\infty f_{2p+1}(x)dx \approx \frac{1-(2p+1)^2}{\pi} \left( -\int_A^\infty \frac{2\sin(2x)}{x}dx - \left[ \frac{\cos 2x}{x} \right]_A^\infty \right) \quad (A18) $$

$$ \approx \frac{1-(2p+1)^2}{\pi} \left( 2\mathrm{Si}(2A) - \pi + \frac{\cos 2A}{A} \right) \quad (A19) $$

• for $g_n$ in 0:

$$ \int_0^\epsilon g_n(x)dx \approx \frac{\epsilon^{2n}}{2n(n!)^2 4^n \sigma} \quad (A20) $$

• for $g_n$ in $+\infty$:

$$ \int_A^\infty g_n(x)dx \approx \frac{4}{\pi} \int_A^\infty \frac{\sin^2\left(x - \frac{n\pi}{2} + \frac{\pi}{4}\right)}{x} e^{-\sigma x} dx \quad (A21) $$

$$ \leq \frac{4}{\pi} \int_A^\infty \frac{1}{x} e^{-\sigma x} dx \quad (A22) $$

$$ \leq \frac{4}{\pi} \frac{e^{-\sigma A}}{\sigma A} \quad (A23) $$

285

So if we take $\sigma A$ sufficiently large (in practise we take $\sigma A \simeq 20$), this part can be neglected.

286

287

The following Table 1 sums up the approximations we used to compute the integrals $I_n - I_1$ and $M_n$.

288

<table><caption>Table 1: Summary of the approximated integrals</caption><thead><tr><th></th><th>f<sub>2p</sub></th><th>f<sub>2p+1</sub></th><th>g<sub>n</sub></th></tr></thead><tbody><tr><td>∫<sub>0</sub><sup>ε</sup></td><td>-ε<sup>2</sup>/8σ</td><td>-ε<sup>2</sup>/8σ</td><td>ε<sup>2n</sup>/(2n(n!)<sup>2</sup>4<sup>n</sup>σ)</td></tr><tr><td>∫<sub>A</sub><sup>∞</sup></td><td>1 - 2/π Si(2A)</td><td>(1-(2p+1)<sup>2</sup>/π) (2Si(2A) - π + cos 2A/A)</td><td>0</td></tr></tbody></table>

289

290

1. Flierl, G.R. On the instability of geostrophic vortices. *Journal of Fluid Mechanics* **1988**, pp. 349–388.

291

292

2. Carton, X.J.; McWilliams, J.C. Nonlinear oscillatory evolution of a baroclinically unstable geostrophic vortex. *Dynamics of Atmospheres and Oceans* **1996**, pp. 207–214.

293

294

3. Baey, J.M.; Carton, X. Vortex multipoles in two-layer rotating shallow-water flows. *Journal of Fluid Mechanics* **2002**, pp. 151–175.

295

296

4. Carton, X. Instability of Surface Quasigeostrophic Vortices. *Journal of the Atmospheric Sciences* **2009**, 66, 1051–1062. doi:10.1175/2008JAS2872.1.

297

298

5. Menesuen, C.; Le Gentil, S.; Marchesiello, P.; Ducousso, N. Destabilization of an Oceanic Meddy-Like Vortex: Energy Transfers and Significance of Numerical Settings. *Journal of Physical Oceanography* **2018**, pp. 1151–1168.

299

300

301

6. Badin, G.; Poulin, F.J. Asymptotic scale-dependent stability of surface quasi-geostrophic vortices: semi-analytic results. *Geophysical & Astrophysical Fluid Dynamics* **2019**, 113, 574–593. doi:10.1080/03091929.2018.1453930.

302

303

304

7. Harvey, B.J.; Ambaum, M.H.P. Perturbed Rankine vortices in surface quasi-geostrophic dynamics. *Geophysical & Astrophysical Fluid Dynamics* **2011**, 105, 377–391. doi:10.1080/03091921003694719.

305

306

307

8. Lapeyre, G. Surface Quasi-Geostrophy. *Fluids* **2017**, p. 28.
308 9. Tulloch, R.; Smith, K.S. Quasigeostrophic Turbulence with Explicit Surface Dynamics: Application to the Atmospheric Energy Spectrum. *Journal of the Atmospheric Sciences* **2009**, *66*, 450–467. doi:10.1175/2008JAS2653.1.

311 10. Smith, K.S.; Bernard, E. Geostrophic turbulence near rapid changes in stratification. *Physics of Fluids* **2013**, *25*, 046601. doi:10.1063/1.4799470.

313 11. Charney, J.G. On the scale of atmospheric motions. *Geofys. Publikasjoner* **1948**, pp. 1–17.

314 12. Held, I.; PierreHumbert, R.T.; Garner, S.T.; Swanson, K.L. Surface quasi-geostrophic dynamics. *Journal of Fluid Mechanics* **1995**, pp. 1–20.

316 13. Gwendal, R.; Lien, H.B.; Patrice, K. Influence of the beta-effect on non-modal baroclinic instability. *Quarterly Journal of the Royal Meteorological Society* **2001**, pp. 1375–1388.

318 14. Fischer, C. Linear Amplification and Error Growth in the 2D Eady Problem with Uniform Potential Vorticity. *Journal of the Atmospheric Sciences* **1998**, *55*, 3363–3380. doi:10.1175/1520-0469(1998)055<3363:LAAEGI>2.0.CO;2.

321 15. Carton, X.; Flierl, G.R.; Perrot, X.; Meunier, T.; Sokolovskiy, M. Explosive instability of geostrophic vortices. Part 1. Baroclinic instability. *Theoretical and Computational Fluid Dynamics* **2010**, pp. 125–130.