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

HAL
open science

# Baroclinic instability in the Eady model for two coupled flows

Armand Vic, Xavier Carton, Jonathan Gula

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

## To cite this version:

Armand Vic, Xavier Carton, Jonathan Gula. Baroclinic instability in the Eady model for two coupled flows. 2024. hal-04826684

HAL Id: hal-04826684

https://univ-paris8.hal.science/hal-04826684v1

Preprint submitted on 9 Dec 2024

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
*Geophysical and Astrophysical Fluid Dynamics*
Vol. 00, No. 00, 00 Month 2023, 1–23

# Baroclinic instability in the Eady model for two coupled flows

Armand VIC†*, Xavier CARTON † and Jonathan GULA †‡

† Univ Brest, CNRS, Ifremer, IRD, Laboratoire d'Océanographie Physique et Spatiale (LOPS), IUEM, 29280, Plouzané, France.

‡ Institut Universitaire de France (IUF), Paris, France

(v4.4 released October 2012)

The linear instability of baroclinic flows in two superimposed and coupled, immiscible fluids is studied theoretically. These flows are westerlies, and they are thermally or mechanically coupled. The fluids are stably stratified, internally and mutually (the upper fluid is lighter than the lower fluid). In each fluid, two-level surface quasi-geostrophy governs the evolution of the perturbed westerly flow. The perturbations are horizontal normal modes. Firstly, the models are not coupled and the flow instability in each fluid is validated separately against the results of the classical Eady model of baroclinic instability. Secondly, the two fluids are thermally and/or mechanically coupled.

With thermal coupling, and for meridionally uniform perturbations, a new mode of instability appears for long waves. This pair of unstable modes converges towards the modes of the uncoupled fluids at medium wavelengths. For perturbations with a non trivial meridional structure, the thermal coupling essentially damps the instability. For an upper flow with a larger deformation radius than in the lower flow, the growth rates of the perturbation are therefore more strongly altered in the former than in the latter. With mechanical coupling, the instability is essentially damped at large to medium scales, while the short-wave cut-off is extended towards smaller waves. When the fluids are both thermally and mechanically coupled, these effects add up. This very idealized study is a first step towards studying more realistic cases.

**Keywords:** surface quasi-geostrophy, two fluid coupling, linearized equations, baroclinic instability, Eady model.

## 1. Introduction

In the Earth atmosphere and oceans, winds and currents carry heat from the Equator to the poles (Richardson 1983, Carton 2010). In the complex atmospheric and oceanic flow patterns, west-east flows (zonal) flows are often energetic (e.g. the Gulf Stream). These strong and narrow currents can become unstable and produce dynamical features with different scales (vortices, filaments). Then the zonal currents interact non-linearly with these features (Fernández-Castro et al. 2020, Vic et al. 2021). For instance, the Antarctic Circumpolar Current (ACC) often becomes unstable and produces mesoscale waves and vortices, a few hundred kilometers in size (Carton 2001). Flow instabilities are one of the processes by which energy is exchanged between scales in the ocean or in the atmosphere (Holland and Haidvogel 1980, Tulloch et al. 2011, Meunier et al. 2013, Duarte et al. 2016).

In the literature of planetary fluid dynamics, the response of oceanic motions to atmospheric forcing, or vice versa (one-way interactions) has been amply studied. Less common are theoretical studies of coupled ocean-atmosphere dynamics (that is, with two-way interactions). Recently, the coupling of oceanic and atmospheric flows has been shown to affect ocean dynamics at the mesoscale or submesoscale (Renault et al. 2016, 2017, 2018). In

<footnote>*Corresponding author. Email: armand.vic@univ-brest.fr</footnote>
particular, these latter studies have shown that oceanic mesoscale motions are weaker with ocean-atmosphere coupling than in the presence of atmospheric forcing only. Therefore, the ocean-atmosphere dynamics with two-way interactions is clearly not the superposition of the two dynamics with one-way interactions; this is due to the nonlinearity of these motions.

Another study (Moulin and Wirth 2014) discovered a new type of barotropic instability occurring from the mechanical coupling of the ocean and the atmosphere. It grows first in the ocean at scales of the baroclinic radius of deformation and then spreads to the atmosphere. In a following study (Moulin and Wirth 2016), the same authors showed that the mechanical coupling of the ocean and the atmosphere leads to kinetic energy absorption by the ocean and conversely, that the ocean variability leaves an imprint in the atmospheric state when the coupling between the two fluids is strong.

The present theoretical study is devoted to understanding how the coupling of two immiscible fluids flowing in the same direction modifies their existing instabilities.

In detail, we investigate how two unstable parallel flows, one in a lighter fluid above and one in a heavier fluid below, evolve linearly when they are coupled. We contrast the coupled flow instability with the uncoupled flow instability. Our westerly flows have a vertical shear of horizontal velocity, necessary for baroclinic instability on a rotating planet. We generalize the Eady model (Eady 1949) to two coupled fluids. The mathematical equations governing the Eady model are derived from the surface quasi-geostrophic dynamics (Bretherton 1966, Held et al. 1995, Lapeyre 2017). To couple our fluids, we propose two different mechanisms. Firstly, we add a thermal relaxation at their interface to mimic a sensible heat flux proportional to the temperature difference between the two fluids; this is called "thermal coupling" hereafter. Secondly, we investigate how a frictional stress between the two fluids influences the linear instabilities. This is called "mechanical coupling" in the following. These formulations of two-fluid coupling, which retain some physical realism, are simple enough to allow an analytical approach. The most unstable perturbations, their intensity and phases, are computed for each fluid. Finally, the two coupling mechanisms are jointly in effect, and the resulting instability is contrasted with that with only one coupling mechanism.

The paper is organized as follows: we present the surface quasi-geostrophic (SQG) equations in section 2. In section 3, we calculate the linearized equations for the instability of thermally coupled westerly flows and we provide an asymptotic analytical solution. Our numerical model of the linear instability problem is validated against the Eady solution for uncoupled fluids. Then, these equations are solved numerically for the thermally coupled model. We conduct a sensitivity study of the growth rates and wavelengths of the perturbations to the physical parameters and the coupling strength. In section 4, we detail the calculations and numerical results for the linear instability of mechanically coupled flows. The results are discussed (section 5) and a conclusion is given (section 6).

## 2. Surface quasi-geostrophic model and equations

### 2.1. Basic equations: the surface quasi-geostrophic model

In this study, we use the surface quasi-geostrophic (SQG) model (Bretherton 1966, Held et al. 1995, Lapeyre 2017). This model applies to stratified fluids, in a rapidly rotating environment, for which the potential vorticity is null in the bulk of the fluid. Potential vorticity is

$$q = [\partial_x^2 + \partial_y^2] \psi + \partial_z \left[ \frac{f_0^2}{N^2} \partial_z \psi \right] \quad (1)$$
where $\psi$ is the streamfunction. $N^2$ is the squared Brunt Väisälä frequency: $N^2 = -(g/\rho_0) \partial\rho/\partial z$.

When potential vorticity is null in the fluid except at its surface and bottom, the flow is governed by surface quasi-geostrophic equations, describing the evolution of $b$, the buoyancy anomaly with respect to a state of rest.

$$ \frac{Db}{Dt} + wN^2 = 0 \quad \text{for } z = 0, H \quad \text{and} \quad b = f_0\partial_z\psi \qquad (2) $$

In this equation $w$ is the vertical velocity at the surface and bottom of the fluid and we have

$$ \frac{Db}{Dt} = \partial_t b + u\partial_x b + v\partial_y b $$

The horizontal velocity is: $u = -\partial_y \psi, v = \partial_x \psi$.

This surface quasi-geostrophic equation is the basis of the Eady model for the baroclinic instability of a single fluid flow. Baroclinic instability affects stratified and vertically sheared flows (Eady 1949).

## 2.2. Basic equations : the four layer case and the various couplings at the interface

In the present study, we consider two superimposed, immiscible fluids, coupled at their interface ($z = 0$) by an exchange of heat (thermal coupling) or by friction (mechanical coupling). Each fluid can undergo baroclinic instability. The global evolution of the coupled fluids will be described by two coupled SQG models. In the upper fluid (that we label "a", for $z > 0$), the streamfunction is $\psi_a(x, y, z, t)$ and the buoyancy anomaly is $b_a = f_0\partial_z\psi_a$. In the lower fluid (that we label "o", for $z < 0$), they are $\psi_o(x, y, z, t)$ and $b_o = f_0\partial_z\psi_o$. We call $H_o$ and $H_a$ the thicknesses of the two fluids. The radii of deformation are the ratios $\sigma_o = \frac{N_a H_o}{f_0}$ in the lower fluid and $\sigma_a = \frac{N_a H_o}{f_0}$ in the upper fluid.

At the upper and lower surfaces of each fluid (levels 1, 2, 3 and 4: top and bottom of the upper and lower fluids, see figure 1), the buoyancy anomalies are:

$$ \psi_a(x, z = H_a, t) = \psi_1(x, y, t), \quad f_0 \frac{\partial \psi_a}{\partial z} \Big|_{z=H_a} = b_1, \qquad (3a) $$

$$ \psi_a(x, z = 0^+, t) = \psi_2(x, y, t), \quad f_0 \frac{\partial \psi_a}{\partial z} \Big|_{z=0} = b_2, \qquad (3b) $$

$$ \psi_o(x, z = 0^-, t) = \psi_3(x, y, t), \quad f_0 \frac{\partial \psi_o}{\partial z} \Big|_{z=0} = b_3, \qquad (3c) $$

$$ \psi_o(x, z = -H_o, t) = \psi_4(x, y, t), \quad f_0 \frac{\partial \psi_o}{\partial z} \Big|_{z=-H_o} = b_4. \qquad (3d) $$

In the absence of forcing, dissipation and of coupling, these anomalies are advected by the local current

$$ \frac{Db_j}{Dt} = \partial_t b_j + u_j \partial_x b_j + v_j \partial_y b_j = 0, \quad j = 1, 4 \qquad (4) $$

When the two fluids are thermally coupled, we assume that the vertical velocity at each
<figure><img src="image_4.png" /><figcaption>Figure 1. Sketch of the 4 layer surface quasi-geostrophic model; the coupling of the upper and lower fluids is achieved between levels 2 and 3. The mean flows at each level are represented.</figcaption></figure>

surface is null; the buoyancies obey the following transport equations at each level:

$$ \frac{Db_1}{Dt} = \frac{Db_4}{Dt} = 0 \qquad (5a) $$

$$ \frac{Db_2}{Dt} = C_o(b_3 - b_2) + F_o \qquad (5b) $$

$$ \frac{Db_3}{Dt} = C_a(b_2 - b_3) + F_a \qquad (5c) $$

where $C_o, C_a$ are the coupling coefficients (determining the sensible heat flux between the two fluids).

When the two fluids are mechanically coupled, friction occurs between levels 2 and 3. Dewar and Flierl (1987) derived the stress-driven coupling between the two fluids. In a stratified quasi-geostrophic model, the influence of the wind on the ocean is achieved via the Ekman pumping (the vertical velocity due to the wind-stress curl). In a surface quasi-geostrophic model, this vertical velocity is easily introduced in the buoyancy equation for the lower fluid as:

$$ \frac{Db_3}{Dt} + w_3 N_o = F_o \qquad (6) $$

where $F_o$ is a supplementary forcing (which renders the mean flow stationary) and $w_3$ is the
Ekman vertical velocity for the lower fluid:

$$w_3 = \frac{1}{\rho_o f_0} [\vec{\nabla} \times \vec{\tau}_2] \cdot \vec{e}_z = \frac{1}{\rho_o f_0} [\partial_x(\tau_2)_y - \partial_y(\tau_2)_x] \quad (7)$$

where $\vec{e}_z$ is the vertical unit vector and $\vec{\tau}_2 = \begin{pmatrix} (\tau_2)_x \\ (\tau_2)_y \end{pmatrix}$ is the stress exerted by the upper fluid flow on the lower one. Generally, the wind stress is a quadratic function of the velocity difference between levels 2 and 3:

$$\vec{\tau}_2 = \rho_a C_d |\vec{u}_2 - \vec{u}_3| (\vec{u}_2 - \vec{u}_3) \quad (8)$$

where $C_d$ is the drag coefficient and $\rho_a$ is the density of the upper fluid. If $|\vec{U}_2| \gg |\vec{U}_3|$, this stress could be linearized as

$$\vec{\tau}_2 = D_a (\vec{u}_2 - \vec{u}_3)$$

where $D_a = \rho_a C_d U_2$.

Note that the stress of the upper on the lower fluid is opposite to that of the lower on the upper fluid. But the vertical velocities divide this stress by the density of the local fluid. Therefore we have $w_2 = -(\rho_o/\rho_a)w_3$. Vertical velocities are much faster in the upper than in the lower fluid, when the former is much lighter.

We assume that the mean flow is parallel and zonal, $V_a = V_o = 0$, as in (Dewar and Flierl 1987). We also assume that the mean flows are steady and thus, that the supplementary forcing exactly balances the stress exerted by each fluid on the other.

In the absence of supplementary forcing, the mean flow would vary with time, and a normal mode analysis would not be appropriate. This case is beyond the scope of the present paper.

## 2.3. Basic state

The basic state depends linearly on $y$ and on $z$ in each fluid:

$$\Psi_a(y, z) = -\Pi_a y - \Lambda_a yz \quad (9a)$$

$$\Psi_o(y, z) = -\Pi_o y - \Lambda_o yz \quad (9b)$$

where $\Lambda_a$, $\Lambda_o$, $\Pi_a$ and $\Pi_o$ are positive constants. The meridional velocities are null $V_a = V_o = 0$ and the zonal velocities increase with $z$:

$$U_a(z) = -\partial_y \Psi_a = \Lambda_a z + \Pi_a \quad (10a)$$

$$U_o(z) = -\partial_y \Psi_o = \Lambda_o z + \Pi_o \quad (10b)$$

or at each level:

$$U_1 = \Pi_a + \Lambda_a H_a, \quad U_2 = \Pi_a, \quad U_3 = \Pi_o, \quad U_4 = \Pi_o - \Lambda_o H_o \quad (11)$$

The mean buoyancies in the four layers are:

$$B_1 = B_2 = -f_0 \Lambda_a y \quad (12a)$$

$$B_3 = B_4 = -f_0 \Lambda_o y \quad (12b)$$

To ensure stationarity of $B_j$, the forcings $F_o$ and $F_a$ are:

$$F_o = -f_0 C_o y (\Lambda_a - \Lambda_o) \quad (13a)$$

$$F_a = f_0 C_a y (\Lambda_a - \Lambda_o) \quad (13b)$$
## 2.4. Perturbations

Adding a normal mode disturbance to the mean flow in each fluid, the total buoyancies are $b_i = B_i + b'_i$ at the four levels $i = 1, 2, 3, 4$, with:

$$b'_i = \beta_i e^{i(kx+ly-\omega t)} \quad (14)$$

From this, we obtain the streamfunction in the two fluids

$$\psi'_j = \varphi_j(z) e^{i(kx+ly-\omega t)} \quad (15)$$

for $j = a, o$ and $\varphi_j$ are functions of $\beta_1, \beta_2, \beta_3, \beta_4$. The assumptions of zero potential vorticity in each fluid yield:

$$-(k^2 + l^2) \varphi_j + \frac{f_0^2}{N_j^2} \partial_z^2 \varphi_j = 0. \quad (16)$$

With $K^2 = k^2 + l^2$, we can write solutions as:

$$\varphi_a(z) = \lambda_a \cosh \left( \frac{N_a}{f_0} K z \right) + \mu_a \cosh \left( \frac{N_a}{f_0} K (z - H_a) \right) \quad (17a)$$

$$\varphi_o(z) = \lambda_o \cosh \left( \frac{N_o}{f_0} K z \right) + \mu_o \cosh \left( \frac{N_o}{f_0} K (z + H_o) \right), \quad (17b)$$

where $\lambda_j, \mu_j$ are computed from conditions (3):

$$\lambda_a = \frac{\beta_1}{s_a}, \quad \mu_a = -\frac{\beta_2}{s_a} \quad (18a)$$

$$\lambda_o = -\frac{\beta_4}{s_o}, \quad \mu_o = \frac{\beta_3}{s_o} \quad (18b)$$

where $s_a = N_a K \sinh \left( \frac{N_a}{f_0} KH_a \right)$, $t_a = N_a K \tanh \left( \frac{N_a}{f_0} KH_a \right)$, $s_o = N_o K \sinh \left( \frac{N_o}{f_0} KH_o \right)$, and $t_o = N_o K \tanh \left( \frac{N_o}{f_0} KH_o \right)$.

# 3. Linear instability of the thermally coupled two-fluid flow

## 3.1. Dispersion relation

We expand the buoyancy equations (5) at first order in the perturbation at the four levels; the left-hand side is:

$$\frac{Db_i}{Dt} = \frac{\partial b'_i}{\partial t} + U_i \frac{\partial b'_i}{\partial x} + v'_i \frac{\partial B_i}{\partial y} \quad (19)$$

and the right-hand side is coupled for the interface layers $i = 1, 2$ and $0$ for the two other layers $i = 1, 4$: the details of the calculation are provided in the first appendix.

The four dispersion relations (A.3), (A.6), (A.9) and (A.12) can be written in a matrix way:

$$M\beta = 0, \text{ where } \beta = \begin{pmatrix} \beta_1 \\ \beta_2 \\ \beta_3 \\ \beta_4 \end{pmatrix} \text{ and}$$
$$ \mathcal{M} = \begin{pmatrix} \omega - k \left( \Pi_a + \Lambda_a \left( H_a - \frac{f_0}{t_a} \right) \right) & \frac{-f_0 k \Lambda_a}{s_a} & 0 & 0 \\ \frac{f_0 k \Lambda_a}{s_a} & \omega - k \left( \Pi_a + \frac{f_0 \Lambda_a}{t_a} \right) + iC_o & -iC_o & 0 \\ 0 & -iC_a & \omega - k \left( \Pi_o - \frac{f_0 \Lambda_o}{t_o} \right) + iC_a & - \frac{f_0 k \Lambda_o}{s_o} \\ 0 & 0 & \frac{f_0 k \Lambda_o}{s_o} & \omega - k \left( \Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right) \right) \end{pmatrix} \quad (20) $$

We want a non trivial solution of the equation $\mathcal{M}\beta = 0$. This means that $\det(\mathcal{M})$ should vanish. Analytically, $\det(\mathcal{M})$ can be computed but we failed to find simple analytical solutions in $\omega$:

$$ \det(\mathcal{M}) = \left[ \omega - k \left( \Pi_a + \Lambda_a \left( H_a - \frac{f_0}{t_a} \right) \right) \right] \left[ \omega - k \left( \Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right) \right) \right] \\ \left[ \omega - k \left( \Pi_a + \frac{f_0 \Lambda_a}{t_a} \right) + iC_o \right] \left[ \omega - k \left( \Pi_o - \frac{f_0 \Lambda_o}{t_o} \right) + iC_a \right] \\ + \frac{f_0^2 k^2 \Lambda_o^2}{s_o^2} \left[ \omega - k \left( \Pi_a + \Lambda_a \left( H_a - \frac{f_0}{t_a} \right) \right) \right] \left[ \omega - k \left( \Pi_a + \frac{f_0 \Lambda_a}{t_a} \right) + iC_o \right] \\ + \frac{f_0^2 k^2 \Lambda_a^2}{s_a^2} \left[ \omega - k \left( \Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right) \right) \right] \left[ \omega - k \left( \Pi_o - \frac{f_0 \Lambda_o}{t_o} \right) + iC_a \right] \\ + C_o C_a \left[ \omega - k \left( \Pi_a + \Lambda_a \left( H_a - \frac{f_0}{t_a} \right) \right) \right] \left[ \omega - k \left( \Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right) \right) \right] \\ + \frac{f_0^4 k^4 \Lambda_a^2 \Lambda_o^2}{s_a^2 s_o^2} $$

$$ \det(\mathcal{M}) = \omega^4 - \omega^3 (\delta_1 + \delta_2 + \delta_3 + \delta_4) \\ + \omega^2 (A + B + C + \delta_1 \delta_2 + \delta_1 \delta_3 + \delta_1 \delta_4 + \delta_2 \delta_3 + \delta_2 \delta_4 + \delta_3 \delta_4) \\ - \omega (A (\delta_1 + \delta_2) + B (\delta_3 + \delta_4) + C (\delta_1 + \delta_4) + \delta_1 \delta_2 \delta_3 + \delta_1 \delta_2 \delta_4 + \delta_1 \delta_3 \delta_4 + \delta_2 \delta_3 \delta_4) \\ + D + A \delta_1 \delta_2 + B \delta_3 \delta_4 + C \delta_1 \delta_4 + \delta_1 \delta_2 \delta_3 \delta_4 $$

where $\delta_1 = k \left( \Pi_a + \Lambda_a \left( H_a - \frac{f_0}{t_a} \right) \right)$, $\delta_2 = k \left( \Pi_a + \frac{f_0 \Lambda_a}{t_a} \right) - iC_o$, $\delta_3 = k \left( \Pi_o - \frac{f_0 \Lambda_o}{t_o} \right) - iC_a$, $\delta_4 = k \left( \Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right) \right)$ are the diagonal coefficient of the following matrix $\mathcal{A}$ defined in (23) and $A = \frac{f_0^2 k^2 \Lambda_o^2}{s_o^2}$, $B = \frac{f_0^2 k^2 \Lambda_a^2}{s_a^2}$, $C = C_o C_a$, $D = \frac{f_0^4 k^4 \Lambda_a^2 \Lambda_o^2}{s_a^2 s_o^2} = AB$.

### 3.2. Long wave $k \to 0$ or short wave $k \gg 1$ limits and finite width perturbations $l \neq 0$

The determinant can easily be computed for long zonal waves $k \to 0$ with finite meridional scale $l \neq 0$: the matrix $\mathcal{M}$ becomes

$$ \lim_{k \to 0} \mathcal{M} = \begin{pmatrix} \omega & 0 & 0 & 0 \\ 0 & \omega + iC_o & -iC_o & 0 \\ 0 & -iC_a & \omega + iC_a & 0 \\ 0 & 0 & 0 & \omega \end{pmatrix} \quad (21) $$

and the only non-real solution in $\omega$ for $\det \lim_{k \to 0} \mathcal{M}$ is:

$$ -i[C_o + C_a] < 0 \quad (22) $$
This solution can easily be obtained by considering a null mean flow (the mean flow terms are multiplied by $ik$), and by subtracting the two linearized equations for the lower level of the atmosphere and for the upper level of the ocean. We obtain $\partial_t(b'_2 - b'_3) = -(C_o + C_a)(b'_2 - b'_3)$. For very long wave perturbations, coupling leads to the damping of the difference between the lower atmospheric waves and the surface oceanic waves, as it will shown later in the numerical computation of the growth rates.

For very short waves $k \to \infty$, with finite $l$, the coupling becomes inefficient (the coupling coefficient is $iC$ in the linearized equations while the other terms are proportional to $k$). We recover two independent Eady problems, one for the atmosphere and one for the ocean, which are stable in the short wave limit. The coupled system is also stable.

### 3.3. Numerical method

Numerically, this problem is solved with eigenvalue routines; indeed $M = \omega \text{Id} - A$. Solutions for $\omega$ such that $\det(M) = 0$ are given by the eigenvalues of $A$ with

$$A = \begin{pmatrix} k \left( \Pi_a + \Lambda_a \left( H_a - \frac{f_0}{t_a} \right) \right) & \frac{f_0 k \Lambda_a}{s_a} & 0 & 0 \\ \frac{-f_0 k \Lambda_a}{s_a} & k \left( \Pi_a + \frac{f_0 \Lambda_a}{t_a} \right) - iC_o & iC_o & 0 \\ 0 & iC_a & k \left( \Pi_o - \frac{f_0 \Lambda_o}{t_o} \right) - iC_a & \frac{f_0 k \Lambda_o}{s_o} \\ 0 & 0 & \frac{-f_0 k \Lambda_o}{s_o} & k \left( \Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right) \right) \end{pmatrix} \quad (23)$$

From the diagonalisation of $A$ we obtain four eigenvalues $\{\omega_1, \omega_2, \omega_3, \omega_4\}$ to which correspond eigenvectors $\{X_1, X_2, X_3, X_4\}$. Each eigenvector has 4 complex coordinates $X_{ij} = A_{ij} \exp(i\phi_{ij})$. Hereafter, we select the two eigenvectors which grow the fastest (having the largest imaginary part of their eigenvalues).

### 3.4. Instability of the uncoupled fluids

When the two fluids are not thermally coupled, we have run our algorithm and we recover the usual curve for baroclinic instability of a two-level SQG flow Eady (1949). The growth rates for two identical fluids, and for an upper fluid with a deformation radius ten times larger than that of the lower fluid, are shown in figure 2. The parameters used here are $l = .0$, $\Pi_a = \Pi_o = 0.0$, $\Lambda_a = \Lambda_o = 1.0$, $R_a = R_0 = 1.0$ for the former case, or $R_1 = 10.0$, $R_o = 1.0$ for the latter case. Obviously, in this latter case, the curve of instability for the upper fluid ("a") is shrunk in wave-numbers, by a factor 10 with respect to that for the lower fluid ("o").

In this uncoupled case, pairs of eigenvectors have complex conjugate eigenvalues, and they are confined in one fluid or in the other. Their vertical phase shift is opposite (in sign) to the vertical shear of mean velocity ($\Lambda$). This is in agreement with the usual Eady model.

Note that $\Pi_a$ and $\Pi_o$ do not influence the instability, and that the growth rates are proportional to $\Lambda_a$ and to $\Lambda_o$.

Note that in the following we will keep $R_o = 1.0$ and $\Lambda_o = 1.0$; this provides a length scale and a time scale $T = 1/\Lambda_o$.

### 3.5. Numerical results on the instability of the coupled system

#### 3.5.1. Influence of the coupling constant at $l = 0$ for two identical fluids

Firstly we consider two identical and coupled fluids with $C_o = C_a = C$. We plot the two largest growth rates with the same parameters as for the uncoupled fluids and we vary $C$. The curve "sig0" corresponds to the uncoupled flow growth rates, while "sig1" and "sig2" show
<figure><img src="image_5.png" /><figcaption>Figure 2. Growth rates of the Eady baroclinic instability for uncoupled fluids; (a) identical fluids, (b) the upper fluid has a deformation radius ten times larger than that of the lower fluid (see all parameters in the text).</figcaption></figure>

the largest and second largest growth rates for coupled flows.

Figure 3 (with $l = 0$) indicates that:

(a) with the coupling, a new mode of instability appears at very long waves; this pair of unstable modes has a weakly and a strongly amplified component. They are not accounted for by the theoretical asymptotic analysis above which relies on $l \neq 0$. No simple asymptotic analysis can be carried out when $l = 0$.

(b) this pair of unstable modes for long waves become more unstable than the modes of the uncoupled system, when the coupling constant increases; this occurs on a growing part of the
<figure><img src="image_6.png" /><figcaption>Figure 3. Growth rates of the Eady baroclinic instability for two identical, thermally coupled fluids with $C_o = C_a$ and $l = 0$. The coupling constant increases are 5 $10^{-4}$ (a) and $10^{-2}$ (b).</figcaption></figure>

unstable wave-number domain.

(c) at $k = 0.1$, and $C = 5 \ 10^{-4}$, the eigenvectors have a nearly constant amplitude vertically (over the two fluids); their vertical phase shift is weak in each fluid. When $C = 10^{-2}$, the phase shift remains small in each fluid, but for the fastest growing mode, a negative phase shift now exists at the interface between the two fluids (on the order of $-\pi/8$).

(d) for medium wavelength perturbations ($k = 1.5$), the previous pair of unstable modes gives way to unstable modes confined in a single fluid (the upper one or the lower one); these latter are the modes of the uncoupled system, weakly damped near the interface, by the coupling.
<figure><img src="image_7.png" /><figcaption>Figure 4. Growth rates of the Eady baroclinic instability for two identical, thermally coupled fluids with $C_o = C_a = 10^{-2}$ and $l = 1$.</figcaption></figure>

Their vertical phase shift is close to $-\pi/4$ in each fluid, and a phase shift of $-\pi/2$ exists between the two fluids (see also Appendix B).

(e) farther in $k$ range, the short-wave cut-off of the instability is hardly modified by the coupling.

Note that, for the rest of this section, $C = 10^{-2}$ corresponds to a weak thermal coupling (associated with a time scale $T_c = 100.0$).

### 3.5.2. Instability of two identical coupled fluids at $l = 1$

Now, we consider the case $l = 1$, again with two identical fluids. Squire's theorem (Squire 1933) indicates that the growth rates in this case should be smaller than those for $l = 0$, which is observed (see below). Figure 4 presents the growth rates for the same flow parameters as for the uncoupled fluid, but with $C = 10^{-2}$. The curves "sig0" correspond to the uncoupled flow growth rates, while "sig1" and "sig2" show the growth rates for coupled flows.

Clearly, the coupling only weakly damps the instability in the two fluids. The numerical results agree with the asymptotic analysis: the system is stable in both the short-wave and the very long wave limits. The coupling only weakly alters the short-wave cut-off of the instability.

### 3.5.3. Linear instability of two different, thermally coupled fluids

First, we consider an upper fluid with a deformation radius ten times larger than that of the lower fluid. We keep $C_o = C_a = C = 10^{-2}$. Figure 5 shows the growth rates of the unstable modes (for the coupled system) compared with those of the uncoupled system. The curves "sigo" and "siga" are the latter, while "sig1" and "sig2" are the former. In this case, the coupling creates four specific modes:

(a) for $k \le 0.15$ a very fast growing mode with growth rates $\sigma \in [0.25, 0.7]$; its phase shift lies between $-\pi/4$ and $-\pi/2$ in the upper fluid; note that a phase shift of $-\pi/2$ occurs at the interface between the two fluids, in relation with the $iC/k$ coupling term.

(b) these two (fast and slowly) growing modes give way to a new pair of unstable modes for $k \in [0.15, 0.25]$ when the growth rates of the fastest growing mode reach those of mode "a".

(c) then, between $k = 0.25$ and $k = 0.6$, the growth rates of these modes join that of the
<figure><img src="image_8.png" /><figcaption>Figure 5. Growth rates of the Eady baroclinic instability for two thermally coupled fluids with $R_o = 1.0$, $R_a = 10.0$, $C_o = C_a = 10^{-2}$ and $l = 0$.</figcaption></figure>

modes of the uncoupled system.

(d) Beyond $k = 0.65$, few differences are observed in the growth rates of this coupled flow, with those of the uncoupled flow, but for a very slight change of the cut-off wave-number.

Second, we consider two fluids with the same deformation radius, but with a vertical shear in the upper fluid tenfold that of the lower fluid (see figure 6. Clearly, the same effect of the thermal coupling is observed: for $k \le 0.5$, the growth rates of the perturbation increase with respect to those of the uncoupled flow, only weakly for the lower fluid, and more strongly for the upper fluid.

# 4. Linear instability of the mechanically coupled ocean-atmosphere flow: analytical and numerical results

## 4.1. Equations for the linear instability of two mechanically coupled SQG fluids

We linearise the coupled buoyancy equations around this mean flow. We have

$$ \partial_t b'_o + U_o \partial_x b'_o + v'_o \partial_y B_o + w'^o_c N_o^2 = 0. $$

The Ekman vertical velocity is, for the perturbation

$$ w'^o_c = \frac{\rho_a C_d}{\rho_o f_0} [U_a - U_o] [\partial_x(v'_a - v'_o) - 2\partial_y(u'_a - u'_o)] \quad (24) $$

Therefore with our normal-mode perturbation, the linearized buoyancy equation with mechanical coupling, is

$$ \partial_t b'_o + U_o \partial_x b'_o + v'_o \partial_y B_o = - \frac{\rho_a C_d N_o^2}{\rho_o f_0} (k^2 + 2l^2) [\psi'_a - \psi'_o] [U_a - U_o] \quad (25) $$

at the ocean surface (layer 3), and

$$ \partial_t b'_o + U_o \partial_x b'_o + v'_o \partial_y B_o = 0 \quad (26) $$
<figure><img src="image_9.png" /><figcaption>Figure 6. Growth rates of the Eady baroclinic instability for two thermally coupled fluids with $R_o = 1.R_a = 1.0$, $\Lambda_a = 10.0$, $\Lambda_o = 1.0$, $C_o = C_a = 10^{-2}$ and $l = 0$.</figcaption></figure>

at the ocean bottom (layer 4), and conversely for the atmosphere, in which the "a" and "o" indices in the formulae above must be permuted.

Now the variables are Fourier decomposed and the equations are divided by $ik$; this provides a pure imaginary coupling of the perturbation equations, as for the thermal coupling. Finally we have to solve the following matrix equation: $\mathcal{M}\beta = 0$ where

$$ \mathcal{M} = \begin{pmatrix} \omega - k \left( \frac{\Pi_a + \Lambda_a (H_a - \frac{f_0}{t_a})}{\frac{f_0 k \Lambda_a}{s_a} - \frac{i M C_o}{s_a}} \right) & \frac{-f_0 k \Lambda_a}{s_a} & 0 & 0 \\ \omega - k \left( \frac{\Pi_a + \frac{f_0 \Lambda_a}{t_a}}{\frac{i M C_o}{t_a}} \right) + \frac{i M C_o}{t_a} & \frac{i M C_o}{t_o} & - \frac{i M C_o}{s_o} & 0 \\ \frac{-i M C_o}{s_a} & \frac{-f_0 k \Lambda_o}{t_o} & \frac{i M C_o}{t_o} & - \frac{f_0 k \Lambda_o}{s_o} \\ 0 & 0 & \omega - k \left( \frac{\Pi_o - \frac{f_0 \Lambda_o}{t_o}}{\frac{f_0 k \Lambda_o}{s_o}} \right) + \frac{i M C_o}{t_o} & \omega - k \left( \frac{\Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right)}{\frac{f_0 k \Lambda_o}{s_o}} \right) \end{pmatrix} \quad (27) $$

where $MC_a = -MC_o = -\frac{\rho_a C_d N_o^2}{\rho_o f_0} (k^2 + 2l^2) [U_2 - U_3]$ are the mechanical coupling coefficients. This is equivalent to find the eigenvalues of the matrix $\mathcal{A}$ where

$$ \mathcal{A} = \begin{pmatrix} k \left( \frac{\Pi_a + \Lambda_a (H_a - \frac{f_0}{t_a})}{\frac{f_0 k \Lambda_a}{s_a} + \frac{i M C_o}{s_a}} \right) & \frac{-f_0 k \Lambda_a}{s_a} & 0 & 0 \\ k \left( \frac{\Pi_a + \frac{f_0 \Lambda_a}{t_a}}{\frac{i M C_o}{t_a}} \right) - \frac{i M C_o}{t_a} & - \frac{i M C_o}{t_o} & \frac{i M C_o}{s_o} & 0 \\ \frac{-i M C_o}{s_a} & \frac{-f_0 k \Lambda_o}{t_o} & - \frac{i M C_o}{t_o} & \frac{f_0 k \Lambda_o}{s_o} + \frac{i M C_o}{s_o} \\ 0 & 0 & k \left( \frac{\Pi_o - \frac{f_0 \Lambda_o}{t_o}}{- \frac{f_0 k \Lambda_o}{s_o}} \right) & k \left( \frac{\Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right)}{\frac{f_0 k \Lambda_o}{s_o}} \right) \end{pmatrix} \quad (28) $$

## 4.2. Asymptotics of growth rates for long zonal waves $k \to 0$ and finite width perturbations $l \neq 0$

The determinant of $\mathcal{M}$ can easily be computed for long zonal waves $k \to 0$ with finite meridional scale $l \neq 0$: the matrix $\mathcal{M}$ becomes
<figure><img src="image_10.png" /><figcaption>Figure 7. Growth rates of the Eady baroclinic instability for two mechanically coupled fluids with $R_o = 1.R_a = 1.0$, $\Lambda_a = \Lambda_o = 1.0$, $\rho_a = \rho_o$, (a) $MC_o = MC_a = 0.01$, (b) $MC_o = MC_a = 0.1$ and $l = 0$.</figcaption></figure>

$$ \lim_{k \to 0} \mathcal{M} = \begin{pmatrix} \omega & 0 & 0 & 0 \\ -\frac{i\widetilde{MC}_o}{\tilde{s}_a} & \omega + \frac{i\widetilde{MC}_o}{\tilde{t}_a} & \frac{i\widetilde{MC}_o}{\tilde{t}_o} & -\frac{i\widetilde{MC}_o}{\tilde{s}_o} \\ -\frac{i\widetilde{MC}_a}{\tilde{s}_a} & \frac{i\widetilde{MC}_a}{\tilde{t}_a} & \omega + \frac{i\widetilde{MC}_a}{\tilde{t}_o} & -\frac{i\widetilde{MC}_a}{\tilde{s}_o} \\ 0 & 0 & 0 & \omega \end{pmatrix} \quad (29) $$

where $\widetilde{MC}_j = \lim_{k \to 0} MC_j = \pm \frac{2\rho_a C_d N_o^2 l^2}{\rho_o f_o} [U_2 - U_3]$, $\tilde{s}_j = \lim_{k \to 0} s_j = N_j l \sinh\left(\frac{N_j l}{f_o} H_j\right)$ and $\tilde{t}_j =$
<figure><img src="image_11.png" /><figcaption>Figure 8. Growth rates of the Eady baroclinic instability for two mechanically coupled fluids with $R_o = 1$, $R_a = 1.0$, $\Lambda_a = \Lambda_o = 1.0$, $\rho_a = \rho_o$, $MC_o = MC_a = 0.1$ and $l = 1$.</figcaption></figure>

$$ \lim_{k \to 0} t_j = N_j l \tanh \left( \frac{N_j}{f_0} l H_j \right) \text{ for } j = \text{a, o. A pure imaginary solution to } \det \left( \lim_{k \to 0} \mathcal{M} \right) \text{ is} \\ -i \left[ \frac{\widetilde{MC}_o}{\tilde{t}_a} + \frac{\widetilde{MC}_a}{\tilde{t}_o} \right] < 0 \qquad (30) $$

So the coupling tends to destabilize the system for long zonal waves. But because $\widetilde{MC}_j \to 0$ as $l \to 0$, the destabilization is weaker as $K \to 0$.

## 4.3. Numerical results for the instability of two mechanically coupled fluids

### 4.3.1. Coupling of two identical fluids

Since $\Pi_a = \Pi_o$ cancels out the coupling, a non-zero $\Pi_a - \Pi_o$ will be used in this section. Firstly we assume that $\Lambda_a = \Lambda_o$, $R_a = R_o$, $\Pi_a - \Pi_o = 1.0$ and $\Delta\rho = \rho_a/\rho_o = 1.0$. We calculate the growth rates of the perturbations for two coupling strengths $MC_o = MC_a = 0.01$ and $MC_o = MC_a = 0.1$.

Contrary to the thermal coupling (C), the mechanical coupling term is not constant; it is $MC/t$ or $MC/s$ where $s$ and $t$ denote the hyperbolic sine or tangent; thus it depends on the wave number. Therefore, the time scale for mechanical coupling depends on the length scale. For $k = 1, l = 0$, $R_a = R_o = 1.0$, the term $(1/s$ or $1/t)$ is on order unity. Then, the coupling time scale is on order of tens to one hundred time units (which is $1/\Lambda_o$).

The mechanical coupling damps the instability in the two fluids and slightly increases the cut-off wave-number. This extension is weak for the weaker coupling but clearly noticeable for the stronger one. The vertical phase shift of the unstable modes, for $MC_o = MC_a = 0.1$, is $-\pi/10$. Each unstable mode is located in one fluid only.

Note that, contrary to the thermal coupling, there is no important change of growth rate for long waves. Indeed, for $l = 0$, the mechanical damping term is proportional to $K^2 = k^2$ which
vanishes more rapidly than the other terms of the linear instability equation (proportional to $k$).

Obviously, this remark does not hold for $l \neq 0$ (see the asymptotic analysis in the subsection above). Figure 8 shows that damping occurs then for long waves.

### 4.3.2. Coupling of two different fluids

Firstly, we modify $\Delta\rho = \rho_a/\rho_o = 0.1$, still assuming that $\Lambda_a = \Lambda_o$, $R_a = R_o$, and $\Pi_a - \Pi_o = 1.0$, with $l = 0$ or $l = 1$. The curves "sig0", "sig1" and "sig2" correspond to the uncoupled flow growth rates, and the largest and second largest growth rates for coupled flows. The instability is more damped in the atmosphere than in the ocean (see figure 9). Again the damping is more efficient for non zero meridional wave-number $l$.

Secondly, we modify the ratio of vertical shears of velocities: we set $\Lambda_a/\Lambda_o = 10.0$, keeping all other parameters equal to those of the identical fluid case. Using $MC_o = MC_a = 0.1$ leads to a 30% damping of the growth rates and a shift of the cut-off wave-number from $k = 2.4$ to $k \sim 3.5$, for $l = 0$. This change is in agreement with those mentioned previously.

Finally, we modify the ratio of deformation radii $R_a/R_o = 10.0$. Figure 10 shows again that the growth rates in the upper and lower fluids are damped by the coupling (a little less so in the upper fluid), with an extension of the short wave cut-off in the lower fluid.

### 4.3.3. Thermal and mechanical coupling of the two fluids

We consider here the case of equal vertical shear in the two fluids $\Lambda_a = \Lambda_o$, but different deformation radii $R_a = 10R_o$ with both thermal $C_a = C_o = 10^{-2}$ and mechanical coupling $MC_a = MC_o = 0.1$. The two fluids have the same density. The meridional wave-number is null $l = 0$. The growth rates (figure 11) are comparable to those of the thermally coupled fluids for long waves, while for shorter waves, the mechanical coupling extends the short wave cut-off.

## 5. Discussion

This study of linear instability of two Eady flows, coupled at their interface, has shown that with thermal coupling, and for zonal perturbations, a pair of unstable modes appears at long waves. This pair of modes converges towards the modes of the uncoupled fluids at medium wavelengths. For perturbations with a non trivial meridional structure, this new mode does not exist and thermal coupling essentially damps the instability, as shown by an asymptotic analysis of the 4th degree equations yielding the growth rates. For an upper flow with a larger deformation radius than in the lower flow, the growth rates of the perturbation are therefore more strongly altered in the former than in the latter. With mechanical coupling, the instability is essentially damped at large to medium scales, while the short-wave cut-off is extended towards smaller waves. When the fluids are both thermally and mechanically coupled, these effects add up. An analytical model of the modal amplitudes and phases (Appendix B) helps interpret the effect of thermal coupling on unstable modes at medium wavelengths. In summary, the analytical models are able to address the small and large $k$ cases, with $l \neq 0$, and the medium $k$ cases for small coupling $C$.

In a very idealized manner, one could apply this study to the atmosphere and ocean, using the following parameters: $R_a = 10.0$, $R_o = 1.0$, $\Delta\rho = \rho_a/\rho_o = 10^{-3}$, $\Lambda_a = 3.0$, $\Lambda_o = 1.0$, $l = 0$. Figure 12 shows the growth rates for the thermally coupled flows, with $C_o = C_a = 10^{-2}$. They bear a notable similarity with those of figure 5. The cases $l = 1$ for thermal coupling and $l = 0$; $l = 1$ for mechanical coupling (not shown here) are similar to those of the previous
<figure><img src="image_12.png" /><figcaption>Figure 9. Growth rates of the Eady baroclinic instability for two mechanically coupled fluids with $R_o = R_a = 1.0$, $\Lambda_a = \Lambda_o = 1.0$, $MC_o = MC_a = 0.1 \rho_a = 0.1 \rho_o$, and (a) $l = 0$, (b) $l = 1$.</figcaption></figure>

subsection.

Compared with previous studies, we do not observe here that coupling leads to a new type of instability, as in (Moulin and Wirth 2014), but only that it modifies an existing one. Also, Pedlosky (1975) found an instability arising from the finite amplitude coupling of the atmosphere and ocean in an idealized quasi-geostrophic model. More work should be carried out to assess the evolution of the perturbation in our system, at finite amplitude.
<figure><img src="image_13.png" /><figcaption>Figure 10. Growth rates of the Eady baroclinic instability for two mechanically coupled fluids with $R_o = 1.0$, $R_a = 10.0$, $\Lambda_a = \Lambda_o = 1.0$, $MC_o = MC_a = 0.1 \rho_a = \rho_o$, and $l = 0$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. Growth rates of the Eady baroclinic instability for two mechanically and thermally coupled fluids with $R_o = 1.0$, $R_a = 10.0$, $\Lambda_a = \Lambda_o = 1.0$, $C_a = C_o = 0.01$, $MC_o = MC_a = 0.1 \rho_a = \rho_o$, and $l = 0$.</figcaption></figure>

Finally, we mention that other types of coupling could have been used: radiative fluxes from each fluid to the other, or latent heat fluxes. Our study being preliminary and idealized, we cannot explore all possible coupling mechanisms.
<figure><img src="image_15.png" /><figcaption>Figure 12. Growth rates of the Eady baroclinic instability for two thermally coupled fluids with $R_o = 1.0, R_a = 10.0, \Lambda_a = 3.0, \Lambda_o = 1.0, C_a = C_o = 0.01, \rho_a = 10^{-3}\rho_o, l = 0$. This case is indicative of parameters for the ocean and the atmosphere.</figcaption></figure>

## 6. Conclusions

As mentioned above, our study is idealized and simple. It could be complemented by numerical nonlinear simulations of the evolution of the coupled flow. With such simulations, more realistic velocity profiles, including in particular a horizontal shear, could be studied. This would broaden our analysis to mixed barotropic-baroclinic instabilities, possibly revealing new effects (as in (Moulin and Wirth 2014, Pedlosky 1975)).

Also the quasi-geostrophic framework does not allow cyclone anticyclone asymmetry. Therefore coupling two two-layer shallow water models would be a natural extension of the present work (and of (Moulin and Wirth 2014)). This would also allow a better comparison with the real coupling of the atmosphere and ocean.

Another situation of interest that can be simulated numerically is the coupling of a parallel atmospheric flow (a zonal wind) with a circular oceanic eddy. The influence of a parallel wind on a circular eddy has been studied in the past, but only as a forcing of the former on the latter (Stern 1965, 1966, Morel and Thomas 2009). The retro-action of the eddy on the wind has also been investigated (Small 2008, Chelton and Xie 2010, Ji et al. 2022). How the two effects interact and modify eddy drift in the ocean, or eddy interaction, still has to be determined. This will be the subject of future work.

## Acknowledgements

The authors thank a referee who suggested to adapt the Davies and Bishop model to the problem at hand. This proved very efficient. Mr Armand Vic wishes to thank Ecole Normale Supérieure de Rennes for a PhD grant. This paper is a contribution to the EUREC4A-OA research program.
## Appendix A: Detailed calculation of the linear instability

The equations for the linear instability are, at the four levels in the thermally coupled model: 1st layer:

$$ -i\omega\beta_1 + (\Pi_a + \Lambda_a H_a)ik\beta_1 - f_0ik\Lambda_a\varphi_a(z = H_a) = 0 \quad (A.1) $$

but

$$ \varphi_a(z = H_a) = \lambda_a \cosh \left( \frac{N_a}{f_0} KH_a \right) + \mu_a = \frac{\beta_1}{t_a} - \frac{\beta_2}{s_a} \quad (A.2) $$

so

$$ \left[ \omega - k \left( \Pi_a + \Lambda_a \left( H_a - \frac{f_0}{t_a} \right) \right) \right] \beta_1 - \frac{f_0 k \Lambda_a}{s_a} \beta_2 = 0 \quad (A.3) $$

2nd layer:

$$ -i\omega\beta_2 + \Pi_a ik\beta_2 - f_0 ik\Lambda_a \varphi_a(z = 0) = C_o (\beta_3 - \beta_2) \quad (A.4) $$

but

$$ \varphi_a(z = 0) = \lambda_a + \mu_a \cosh \left( \frac{N_a}{f_0} KH_a \right) = \frac{\beta_1}{s_a} - \frac{\beta_2}{t_a} \quad (A.5) $$

so

$$ \frac{f_0 k \Lambda_a}{s_a} \beta_1 + \left[ \omega - k \left( \Pi_a + \frac{f_0 \Lambda_a}{t_a} \right) + iC_o \right] \beta_2 - iC_o \beta_3 = 0 \quad (A.6) $$

3rd layer:

$$ -i\omega\beta_3 + \Pi_o ik\beta_3 - f_0 ik\Lambda_o \varphi_o(z = 0) = C_o (\beta_2 - \beta_3) \quad (A.7) $$

but

$$ \varphi_o(z = 0) = \lambda_o + \mu_o \cosh \left( \frac{N_o}{f_0} KH_o \right) = \frac{\beta_3}{t_o} - \frac{\beta_4}{s_o} \quad (A.8) $$

so

$$ -iC_o \beta_2 + \left[ \omega - k \left( \Pi_o - \frac{\Lambda_o f_0}{t_o} \right) + iC_o \right] \beta_3 - \frac{k f_0 \Lambda_o}{s_o} \beta_4 = 0 \quad (A.9) $$

4th layer:

$$ -i\omega\beta_4 + (\Pi_o - \Lambda_o H_o) ik\beta_4 - f_0 ik\Lambda_o \varphi_o(z = -H_o) = 0 \quad (A.10) $$

but

$$ \varphi_o(z = -H_o) = \lambda_o \cosh \left( \frac{N_o}{f_0} KH_o \right) + \mu_o = \frac{\beta_3}{s_o} - \frac{\beta_4}{t_o} \quad (A.11) $$

so

$$ \frac{f_0 k \Lambda_o}{s_o} \beta_3 + \left[ \omega - k \left( \Pi_o + \Lambda_o \left( \frac{f_0}{t_o} - H_o \right) \right) \right] \beta_4 = 0 \quad (A.12) $$

## Appendix B: An analytical model for the amplitudes and phases of the unstable modes in the thermally coupled case

### B.1. Equations

Here we follow the nonlinear model of coupled ODE's for the amplitudes and phases of unstable waves in the Eady problem, devised by Davies and Bishop (1994). We will only remain in the
context of normal modes, and we delay the study of singular modes to further work.

We assume two identical fluids with $\Lambda_o = \Lambda_a = \Lambda$, $H_a = H_o = H$, $N_a = N_o = N$, $\mu = NK/f_0$ where K is the isotropic wavenumber, and $\gamma = 1/\sinh(\mu H)$). With our notations, the buoyancy perturbation at each level is

$$b'_{j}(x, y, t) = A_{j}(t) \cos(kx + \phi_{j}(t)) \exp(ily)$$

where $j = 1, 2, 3, 4$ s the level index. In the following, we will provide evolution equations for the $A_j$'s and the $\phi_j$'s.

The perturbation streamfunction is

$$\psi'_{1}(x, y, t) = \frac{\gamma}{\mu} [A_{1}(t) \cosh(\mu H) \cos(kx + \phi_{1}) - A_{2} \cos(kx + \phi_{2})] \exp(ily)$$

$$\psi'_{2}(x, y, t) = \frac{\gamma}{\mu} [A_{1}(t) \cos(kx + \phi_{1}) - A_{2} \cosh(\mu H) \cos(kx + \phi_{2})] \exp(ily)$$

$$\psi'_{3}(x, y, t) = \frac{\gamma}{\mu} [A_{3}(t) \cosh(\mu H) \cos(kx + \phi_{3}) - A_{4} \cos(kx + \phi_{4})] \exp(ily)$$

$$\psi'_{4}(x, y, t) = \frac{\gamma}{\mu} [A_{3}(t) \cos(kx + \phi_{3}) - A_{4} \cosh(\mu H) \cos(kx + \phi_{4})] \exp(ily)$$

By adequately multiplying these equations by $\sin(kx + \phi_j)$ or by $\cos(kx + \phi_j)$ and integrating over a wavelength, we obtain the following set of coupled, nonlinear ODE's

$$dA_{1}/dt = -\frac{\Lambda\gamma k}{\mu}A_{2} \sin(\phi_{1} - \phi_{2})$$

$$d\phi_{1}/dt = -k\Lambda H + \frac{\Lambda\gamma k}{\mu} [\cosh(\mu H) - \frac{A_{2}}{A_{1}}\cos(\phi_{1} - \phi_{2})]$$

$$dA_{2}/dt = -\frac{\Lambda\gamma k}{\mu}A_{1} \sin(\phi_{1} - \phi_{2}) - C [A_{2} - A_{3} \cos(\phi_{2} - \phi_{3})]$$

$$d\phi_{2}/dt = \frac{\Lambda\gamma k}{\mu} [-\cosh(\mu H) + \frac{A_{1}}{A_{2}}\cos(\phi_{1} - \phi_{2})] - C \frac{A_{3}}{A_{2}} \sin(\phi_{2} - \phi_{3})$$

$$dA_{3}/dt = -\frac{\Lambda\gamma k}{\mu}A_{4} \sin(\phi_{4} - \phi_{3}) - C [A_{3} - A_{2} \cos(\phi_{2} - \phi_{3})]$$

$$d\phi_{3}/dt = \frac{\Lambda\gamma k}{\mu} [\cosh(\mu H) - \frac{A_{4}}{A_{3}}\cos(\phi_{3} - \phi_{4})] - C \frac{A_{2}}{A_{3}} \sin(\phi_{2} - \phi_{3})$$

$$dA_{4}/dt = -\frac{\Lambda\gamma k}{\mu}A_{3} \sin(\phi_{3} - \phi_{4})$$

$$d\phi_{4}/dt = k\Lambda H + \frac{\Lambda\gamma k}{\mu} [-\cosh(\mu H) + \frac{A_{3}}{A_{4}}\cos(\phi_{3} - \phi_{4})]$$

## B.2. Solutions

The objective now is to find approximate solutions to this system of 8 equations by assuming that the ratio of the coupling constant to the growth rate is weak (on the order of $\epsilon \ll 1$). To do so, we investigate two cases of the phase shift at the interface between the two fluids,
$\phi_2 - \phi_3 = \pi/2$ and $\phi_2 - \phi_3 = 0$.

Start with the first case, one can show that equations for levels (1,2) and (3,4) decouple, and that they obey the same system of equations, Hence we can assume that $A_2 = A_3 = A_1 = A_4$ by symmetry. Writing $\sigma = \Lambda\gamma/\mu$, $\phi = \phi_1 - \phi_2$ and $2\hat{\phi} = \phi_1 + \phi_2$, we have first the equation

$$d^2 A_2/dt^2 - C dA_2/dt - (\sigma k \sin(\phi))^2 A_2 = 0$$

We set $\chi = \sigma k \sin(\phi)$. With $A_2(t) = A_0 \exp(at)$, we obtain the equation $a^2 - Ca - \chi^2 = 0$. Assuming a weak coupling, we set $C = 2\chi\epsilon$. A short expansion shows that $a = \chi(1 - \epsilon)$. This indicates that the growth rate of the unstable mode is weakly decreased by the coupling.

Then we obtain that $A_2/A_1 = 1 - \epsilon$, that is, the coupling weakly damps the perturbation near the interface. Note that these asymptotic results cannot hold if $k$ is small because $C$ must be small in front of $\chi$.

More information can be gained by analyzing the equations for $\phi$ and for $\hat{\phi}$. After removing the part independent of the coupling, and calling $\delta\phi$ and $\delta\hat{\phi}$ the phase differences, we obtain

$$d(\delta\hat{\phi})/dt = -\epsilon(\chi + \sigma k \cos(\phi))$$

or, if $A_3 = 0$,

$$d(\delta\hat{\phi})/dt = -\epsilon\sigma k \cos(\phi)$$

that is, the additive phase, or propagative phase, is slowed down by the coupling. We also obtain

$$d(\delta\phi)/dt = -2\epsilon\chi$$

or, if $A_3 = 0$,

$$d(\delta\phi)/dt = 0$$

showing that the subtractive phase, or phase difference, is also reduced by the coupling. Thus, the instability is lessened.

In summary, and as observed in the numerical calculations, the effect of the coupling is to damp the growth rates, the amplitude of the wave at the interface, to slow down the waves and to slightly reduce their vertical phase shift, for finite $k$ and in the case of weak coupling. No simple conclusion can be drawn in more complex situations.

Now, we end up with the second case. The equations are (assuming again identical upper and lower fluids)

$$dA_1/dt = -\sigma k A_2 \sin(\phi); \quad dA_2/dt = -\sigma k A_1 \sin(\phi) - C(A_2 - A_3)$$

$$dA_3/dt = -\sigma k A_4 \sin(\phi) - C(A_3 - A_2); \quad dA_4/dt = -\sigma k A_3 \sin(\phi).$$

We can only infer from this, that in the general case the sum of the perturbation amplitudes obeys the same equation as in the uncoupled case, and more particularly, that the coupling will vanish if the unstable waves have the same amplitude at levels 2 and 3.

References

Bretherton, F.P., Critical layer instability in baroclinic flows. *Q. J. R. Meteorol. Soc.*, 1966, **92**, 325–334.

Carton, X.J., Hydrodynamical modelling of oceanic vortices. *Surv. Geophys.*, 2001, **22**, 179–263.

Carton, X.J., Oceanic Vortices; in *Fronts, Waves and Vortices in Geophysical Flows*, edited by J.B. Flor, Vol. Lecture Notes in Physics, 805, 2010, pp. 61–108.

Chelton, D.B. and Xie, S.P., Coupled ocean-atmosphere interaction at oceanic mesoscales. *Oceanography*, 2010, **23**, 52–69.
Davies, H. and Bishop, C., Eady edge waves and rapid development. *J. Atmos. Sci.*, 1994, **51**, 1930—1946.

Dewar, W.K. and Flierl, G.R., Some Effects of the Wind on Rings. *J. Phys. Oceanogr.*, 1987, **17**, 1653—1667.

Duarte, R.M., Carton, X.J. and Poulin, F.J., The dynamics of a meandering coastal jet in the lee of a cape. *Reg. Chaot. Dyn.*, 2016, **21**, 274–290.

Eady, E.T., Long Waves and Cyclone Waves. *Tellus*, 1949, **1**, 33–52.

Fernández-Castro, B., Evans, D.G., Frajka-Williams, E., Vic, C. and Naveira-Garabato, A.C., Breaking of Internal Waves and Turbulent Dissipation in an Anticyclonic Mode Water Eddy. *J. Phys. Oceanogr.*, 2020, **50**, 1893–1914.

Held, I., PierreHumbert, R.T., Garner, S.T. and Swanson, K.L., Surface quasi-geostrophic dynamics. *J. Fluid Mech.*, 1995, pp. 1–20.

Holland, W.R. and Haidvogel, D.B., A parameter study of the mixed instability of idealized ocean currents. *Dyn. Atmos. Oceans*, 1980, **4**, 185–215.

Ji, J., Dong, C., Liu, X., Liu, T., Yu, Y., Lim Kam Sian, K.T., Bethel, B.J. and Zhao, H., Influence of oceanic mesoscale eddy on the atmospheric boundary layer based on an idealized model. *Deep-Sea Res. Part II*, 2022, **202**, 105146.

Lapeyre, G., Surface Quasi-Geostrophy. *Fluids*, 2017, **2**, 1–28.

Meunier, T., Carton, X.J. and Duarte, R., Influence of a deep flow on a surface boundary current. *Geophys. Astrophys. Fluid Dyn.*, 2013, **107**, 277–303.

Morel, Y. and Thomas, L.N., Ekman drift and vortical structures. *Ocean Model.*, 2009, **27**, 185–197.

Moulin, A. and Wirth, A., A Drag-Induced Barotropic Instability in Air-Sea Interaction. *J. Phys. Oceanogr.*, 2014, pp. 733–741.

Moulin, A. and Wirth, A., Momentum Transfer Between an Atmospheric and an Oceanic Layer at the Synoptic and the Mesoscale: An Idealized Numerical Study. *Bound.-Layer Meteorol.*, 2016, **3**, 551–568.

Pedlosky, J., The development of thermal anomalies in a coupled ocean-atmosphere model. *J. Atmos. Sci.*, 1975, **32**, 1501—1514.

Renault, L., Marchesiello, P., Masson, S. and McWilliams, J.C., Remarkable Control of Western Boundary Currents by Eddy Killing, a Mechanical Air-Sea Coupling Process. *Geophys. Res. Lett.*, 2017, **46**, 2743–2751.

Renault, L., McWilliams, J.C. and Gula, J., Dampening of Submesoscale Currents by Air-Sea Stress Coupling in the Californian Upwelling System. *Sci. Rep.*, 2018, **8**, 1–7.

Renault, L., Molemaker, J.M., McWilliams, J.C., Shchepetkin, A.F., Lemarié, F., Chelton, D., Illig, S. and Hall, A., Modulation of Wind Work by Oceanic Current Interaction with the Atmosphere. *J. Phys. Oceanogr.*, 2016, **46**, 1685–1704.

Richardson, P.L., Gulf-Stream rings; in *Eddies in Marine Science*, edited by A.R. Robinson, Vol. Topics in Atmospheric and Oceanographic Sciences, 1983, pp. 19–45.

Small, R.J., Air–sea interaction over ocean fronts and eddies. *Dyn. Atmos. Oceans*, 2008, **45**, 274–319.

Squire, H.B., On the stability for three-dimensional disturbances of viscous fluid flow between parallel walls. *Proc. R. Soc. Lond. A*, 1933, **142**, 621–628.

Stern, M.E., Interaction of a uniform wind stress with a geostrophic vortex. *Deep-Sea Res.*, 1965, **12**, 355–367.

Stern, M.E., Interaction of a uniform wind stress with hydrostatic eddies. *Deep-Sea Res.*, 1966, **13**, 193–203.

Tulloch, R., Marshall, J., Hill, C. and Smith, K.S., Scales, Growth Rates, and Spectral Fluxes of Baroclinic Instability in the Ocean. *J. Phys. Oceanogr.*, 2011, **41**, 1057–1076.

Vic, A., Carton, X. and Gula, J., The Interaction of Two Unsteady Point Vortex Sources in a Deformation Field in 2D Incompressible Flows. *Reg. Chaot. Dyn.*, 2021, **26**, 618–646.