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

HAL
open science

# A unified formulation of quasi-geostrophic and shallow water equations via projection

Louis Thiry, Long Li, Etienne Mémin, Guillaume Roullet

► To cite this version:

Louis Thiry, Long Li, Etienne Mémin, Guillaume Roullet. A unified formulation of quasi-geostrophic and shallow water equations via projection. Journal of Advances in Modeling Earth Systems, 2024, 10.1029/2024MS004510. http://hal-04673756v2

HAL Id: hal-04673756

https://hal.science/hal-04673756v2

Submitted on 27 Sep 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A unified formulation of quasi-geostrophic and shallow water equations via projection

Louis Thiry<fnref n="1" />, Long Li<fnref n="2" />, Etienne Mémin<fnref n="2" />, and Guillaume Roullet<fnref n="3" />

## Key Points:

* Introduces QG equations as a projection of RSW equations using the same prognostic variables.

* Allows consistent nested numerical methods for both QG and RSW equations.

* Successfully tested on vortex shear instability, double-gyre circulation, and a simplified North Atlantic configuration.

<footnote n="1">¹ANGE Team, INRIA Paris, Laboratoire Jacques-Louis Lions, Sorbonne University, France</footnote>

<footnote n="2">²ODYSSEY Team, INRIA Rennes, IRMAR, University of Rennes I, France</footnote>

<footnote n="3">³Laboratoire d'Océanographie Physique et Spatiale, Université de Bretagne Occidentale, Brest, France</footnote>
**Abstract**

This paper introduces a unified model for layered rotating shallow-water (RSW) and quasi-geostrophic (QG) equations, based on the intrinsic relationship between these two sets of equations. We propose a novel formulation of the QG equations as a projection of the RSW equations. This formulation uses the same prognostic variables as RSW, namely velocity and layer thickness, thereby restoring the proximity of these two sets of equations. It provides direct access to the ageostrophic velocities embedded within the geostrophic velocities resolved by the QG equations. This approach facilitates the study of differences between QG and RSW using a consistent numerical discretization. We demonstrate the effectiveness of this formulation through examples including vortex shear instability, double-gyre circulation, and a simplified North Atlantic configuration.

**Plain Language Summary**

In this paper, we present a straightforward way to connect two important sets of ocean equations: the layered rotating shallow-water (RSW) equations and the quasi-geostrophic (QG) equations. We consider a unified method to formulate the QG equations as a projection of the RSW equations. This method uses the same variables as the RSW equations, making it easier to understand how these two sets of equations relate to each other. Our approach also provides direct access to important velocity information that is usually hidden in the QG equations. This allows us to study the differences between the QG and RSW equations using the same numerical techniques. We show the benefits of our model with two examples: a vortex shear instability and a double-gyre configuration. This work provides a useful tool for the understanding of ocean dynamics. Developing a unified numerical framework for nested models can greatly simplify oceanographic modeling and enhance accuracy.

**1 Introduction**

Large-scale ocean models offer a natural trade-off between complexity and realism. On one side, the primitive equations (PE) are sufficiently realistic for climate simulations and real-world data assimilation. However, they are quite complex as they describe the coupling between the equations of motion (mass and momentum) and the conservation of tracers via the equation of state following thermodynamic laws. On the other side, the barotropic planetary geostrophic equations can be efficiently solved with a single prognostic variable (the thickness), yet they ignore thermodynamics and vertical variations.

Between these extremes, there exist approximate models of the PE, such as the multi-layer rotating shallow-water (RSW) and the multi-layer quasi-geostrophic (QG) models. The latter can be either derived from the former using an asymptotic approach (Pedlosky, 2013; Vallis, 2017) or considered as a vertical discretization of the continuously stratified QG model derived from the PE. These layered models describe the dynamics of vertically stratified flow in isentropic (or isopycnal) coordinates and only require solving the horizontal momentum and mass equations. For instance, the multi-layer RSW model can fairly reproduce the ocean dynamics in the Gulf Stream region with only five vertical levels (Hurlburt & Hogan, 2000). The multi-layer QG equations are widely adopted for the development of deterministic mesoscale eddy parameterizations (Marshall et al., 2012; Jansen & Held, 2014; Bachman et al., 2017; Uchida et al., 2022) as well as stochastic ones (Berloff, 2005; Grooms et al., 2015; Zanna et al., 2017; Bauer et al., 2020; Li et al., 2023). These two sets of equations can be used as research tools since their numerical integration is light enough to run on laptop computers.

Although the multi-layer QG equations are formally derived from the multi-layer RSW equations, their relationship is still not clear from a numerical point of view. Indeed, they are usually written with different prognostic variables. The prognostic variables of the RSW system are the horizontal velocity (**u**, **v**) and the layer thickness **h**, whereas the QG system
is typically formulated with potential vorticity $\mathbf{q}$ as the single prognostic variable, from which one can diagnose the streamfunction $\psi$ and the pressure $\mathbf{p}$. This difference breaks the conceptual continuity between these two equation sets in the ocean model hierarchy and brings some undesirable consequences in practice. For instance, many eddy parameterizations (e.g., Bachman, 2019; Li et al., 2023) yield different formulations and/or discretizations when applied to the RSW model using the horizontal velocity ($\mathbf{u}, \mathbf{v}, \mathbf{h}$) and layer thickness as the prognostic variables, or to the QG model using potential vorticity $\mathbf{q}$ as the prognostic variable. Moreover, from a practical point of view, it is not straightforward to use the same discrete schemes for these two different models in order to compare them under the same configuration.

In this work, we propose to reformulate the multi-layer QG model using horizontal velocity and layer thickness ($\mathbf{u}, \mathbf{v}, \mathbf{h}$) as prognostic variables. This unified reformulation is achieved by expressing the QG equations as a projection of the RSW equations. We also present a numerical algorithm to integrate the QG equations corresponding to the proposed formulation. This projection approach allows for the construction of a QG discretization on top of any RSW discretization. Moreover, this formulation provides direct access to the ageostrophic velocity, which is hidden in the standard QG model formulation but contributes to QG dynamics.

Similarly, the QG dynamics can be understood as a projection in the space of normal-mode variables, as already proposed by Leith (1980); Salmon (1998); Saujani and Shepherd (2006). This projection reveals that the slow modes correspond to the linearized potential vorticity of the RSW model, while the fast modes are associated with the ageostrophic curl and divergence components. However, this formal analysis relies on Fourier transformation, which assumes periodic solutions and thereby limits numerical development. In our formulation, the QG equations are considered as a projection of the RSW model in physical space. This approach is independent of the model configuration, such as the shape of the boundary. It provides a consistent framework for developing the corresponding numerical scheme in realistic basins.

For the numerical experiments, we adopt the RSW discretization proposed by Roullet and Gaillard (2022), which relies on the vector invariant formulation and an advanced high-order WENO advection scheme. The discretized QG model uses the same dynamical core as the multi-layer RSW model, with the only modification being the addition of the projection operator. We developed a compact, efficient, and CPU/GPU portable Python code using the PyTorch library (Paszke et al., 2019). We first use a simple test case of vortex shear instability to investigate the similarities and differences of the QG and RSW solutions according to different Rossby numbers. We then study the solutions produced by QG and RSW in an idealized wind-driven double-gyre configuration. This final configuration is also implemented in the North Atlantic basin with simplified lateral and vertical boundary conditions.

This paper is organized as follows. In Section 2, we briefly recall the derivation of the multi-layer RSW and QG equations. In Section 3, we present our projected QG formulation. In Section 4, we numerically test our formulation on different configurations. We conclude and discuss further perspectives in Section 5.

## 2 Multi-layer RSW and QG equations

In this section, we first review the governing equations of the RSW system. We then briefly explain the QG scaling of the RSW equations and present the standard multi-layer QG equations, using potential vorticity $\mathbf{q}$ as the prognostic variable.
## 2.1 Multi-layer RSW equations

The stratification of a multi-layer RSW model consists of a stack of n isopycnal layers, as illustrated in Figure 1 with n = 3 layers. By convention, we index the layers from i = 1 for the top layer to i = n for the bottom layer. These layers have a uniform reference thickness $H_i$ and a density $\rho_i$. As shown in Figure 1, the total thickness of a layer i is the sum of the reference thickness $H_i$ and the thickness anomaly $h_i(x, y)$. The vertical interface displacement $\eta_i(x, y)$ and the hydrostatic pressure $p_i(x, y)$ are given by

$$ \eta_i = \sum_{j=i}^{n} h_j, \quad p_i = \rho_1 \sum_{j=1}^{i} g'_j \eta_j. \qquad (1) $$

For $i > 1$, the reduced gravities are $g'_i = g(\rho_i - \rho_{i-1})/\rho_1$, and for the top layer $g'_1 = g$.

Using the following vector notation

$$ \mathbf{H} = (H_1, \dots, H_n), \qquad (2a) $$

$$ \mathbf{h} = (h_1(x, y), \dots, h_n(x, y)), \qquad (2b) $$

$$ \boldsymbol{\eta} = (\eta_1(x, y), \dots, \eta_n(x, y)), \qquad (2c) $$

$$ \mathbf{p} = (p_1(x, y), \dots, p_n(x, y)), \qquad (2d) $$

we can rewrite the linear relations between $\mathbf{h}, \eta$ and $\mathbf{p}$ in a compact form

$$ \mathbf{h} = \mathrm{diag}(\mathbf{H})A\mathbf{p}, \qquad (3a) $$

$$ \eta = C\mathbf{h}, \qquad (3b) $$

$$ \mathbf{p} = M\mathbf{h}, \qquad (3c) $$

with the matrices

$$ A = \frac{1}{\rho_1} \begin{bmatrix} \frac{1}{H_1 g'_1} + \frac{1}{H_1 g'_2} & \frac{-1}{H_1 g'_2} & \cdot & \cdot \\ \frac{-1}{H_2 g'_2} & \frac{1}{H_2} \left( \frac{1}{g'_2} + \frac{1}{g'_3} \right) & \frac{-1}{H_2 g'_3} & \cdot \\ \cdot & \cdot & \cdot & \cdot \\ \cdot & \cdot & \frac{-1}{H_n g'_n} & \frac{1}{H_n g'_n} \end{bmatrix}, \qquad (3d) $$

$$ B = \rho_1 \begin{bmatrix} g'_1 & 0 & 0 & \dots & 0 \\ g'_1 & g'_2 & 0 & \dots & 0 \\ \cdot & \cdot & \cdot & \dots & \cdot \\ g'_1 & g'_2 & g'_3 & \dots & g'_n \end{bmatrix}, \quad C = \begin{bmatrix} 1 & 1 & \dots & 1 & 1 \\ 0 & 1 & \dots & 1 & 1 \\ \cdot & \cdot & \cdot & \dots & \cdot \\ 0 & 0 & \dots & 0 & 1 \end{bmatrix}, \qquad (3e) $$

and $M = BC$. We then introduce the velocity components ($\mathbf{u}, \mathbf{v}$), the kinetic energy $\mathbf{k}$, and the relative vorticity $\omega$:

$$ \mathbf{u} = (u_1(x, y), \dots, u_n(x, y)), \qquad (4a) $$

$$ \mathbf{v} = (v_1(x, y), \dots, v_n(x, y)), \qquad (4b) $$

$$ \mathbf{k} = (\mathbf{u}^2 + \mathbf{v}^2)/2, \qquad (4c) $$

$$ \omega = \partial_x \mathbf{v} - \partial_y \mathbf{u}. \qquad (4d) $$

With these variables, the multi-layer RSW equations read

$$ \partial_t \mathbf{u} = (\omega + f)\mathbf{v} - \partial_x(\mathbf{p} + \mathbf{k}), \qquad (5a) $$

$$ \partial_t \mathbf{v} = -(\omega + f)\mathbf{u} - \partial_y(\mathbf{p} + \mathbf{k}), \qquad (5b) $$

$$ \partial_t \mathbf{h} = -\mathbf{H}(\partial_x \mathbf{u} + \partial_y \mathbf{v}) - \partial_x(\mathbf{u}\mathbf{h}) - \partial_y(\mathbf{v}\mathbf{h}), \qquad (5c) $$

where $f$ is the Coriolis parameter. These equations can be formulated in the compact form

$$ \partial_t \mathbf{X} = F(\mathbf{X}), \qquad (6a) $$
with $\mathbf{X} = (\mathbf{u}, \mathbf{v}, \mathbf{h})^T$, the state variable, and

$$F \begin{pmatrix} \mathbf{u} \\ \mathbf{v} \\ \mathbf{h} \end{pmatrix} = \begin{pmatrix} (\partial_x \mathbf{v} - \partial_y \mathbf{u} + f)\mathbf{v} - \partial_x (M\mathbf{h} + (\mathbf{u}^2 + \mathbf{v}^2)/2) \\ -(\partial_x \mathbf{v} - \partial_y \mathbf{u} + f)\mathbf{u} - \partial_y (M\mathbf{h} + (\mathbf{u}^2 + \mathbf{v}^2)/2) \\ -\mathbf{H}(\partial_x \mathbf{u} + \partial_y \mathbf{v}) - \partial_x (\mathbf{u}\mathbf{h}) - \partial_y (\mathbf{v}\mathbf{h}) \end{pmatrix}, \quad (6b)$$

the RSW model operator.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Vertical cross-section of a three-layer RSW model displaying layer thicknesses $H_i + h_i$ and interface heights $\eta_i$</figcaption></figure>

## 2.2 QG scaling of RSW equations

The QG model is derived from the RSW model under two scaling assumptions: $Ro \ll 1$ and $Bu \sim 1$ (J. McWilliams, 2006; Zeitlin, 2018), where $Ro$ is the Rossby number and $Bu$ is the Burger number (ratio of the deformation Radius $L_d = \sqrt{gH}/f$ to the scale of motion $L$). The consistency of the QG scaling necessitates a beta-plane approximation for the Coriolis parameter $f = f_0 + \beta y$, where $\beta$ is the meridional gradient of the Coriolis parameter. The condition $Ro \ll 1$ implies that the velocity is near the geostrophic balance, specifically that the ageostrophic correction is $O(Ro)$. We define $(\mathbf{u}_g, \mathbf{v}_g)$ as the geostrophic velocity that adheres the geostrophic balance

$$-f_0 \mathbf{u}_g = \partial_y \mathbf{p}, \quad f_0 \mathbf{v}_g = \partial_x \mathbf{p}, \qquad (7a)$$

and decompose the velocity into the geostrophic and the ageostrophic components $(\mathbf{u}_a, \mathbf{v}_a)$:

$$\mathbf{u} = \mathbf{u}_g + \mathbf{u}_a, \quad \mathbf{v} = \mathbf{v}_g + \mathbf{v}_a. \qquad (7b)$$

The conditions $Bu \sim 1$ and $Ro \ll 1$ imply that $\mathbf{h}/\mathbf{H} \sim O(Ro)$. Following the approach of Holton (1973); J. McWilliams (2006); Cushman-Roisin and Beckers (2011); Zeitlin (2018), we do not decompose $\mathbf{p}, \mathbf{h}, \eta$ into a geostrophic and an ageostrophic parts, as this would require an additional assumption to define the ageostrophic component of the stratification, which is unnecessary for deriving the QG model. Thus, there is a single mass variable and two velocity components per layer.

Under these assumptions, the multi-layer QG equations can be expressed in terms of three equations for the prognostic variables $(\mathbf{u}_g, \mathbf{v}_g, \mathbf{h})$ with the right-hand side terms that depend on the ageostrophic velocities $(\mathbf{u}_a, \mathbf{v}_a)$:

$$\partial_t \mathbf{u}_g = (\omega_g + f_0 + \beta y)\mathbf{v}_g + f_0 \mathbf{v}_a - \partial_x (M\mathbf{h} + \mathbf{k}_g), \quad (8a)$$

$$\partial_t \mathbf{v}_g = -(\omega_g + f_0 + \beta y)\mathbf{u}_g - f_0 \mathbf{u}_a - \partial_y (M\mathbf{h} + \mathbf{k}_g), \quad (8b)$$

$$\partial_t \mathbf{h} = -\mathbf{H}(\partial_x \mathbf{u}_g + \partial_y \mathbf{v}_g) - \partial_x (\mathbf{u}_g \mathbf{h}) - \partial_y (\mathbf{v}_g \mathbf{h}) - \mathbf{H}(\partial_x \mathbf{u}_a + \partial_y \mathbf{v}_a). \quad (8c)$$
Note that using the geostrophic balance (7a), one can simplify these equations to

$$ \partial_t \mathbf{u}_g = (\omega_g + \beta y)\mathbf{v}_g + f_0 \mathbf{v}_a - \partial_x \mathbf{k}_g, \quad (9a) $$

$$ \partial_t \mathbf{v}_g = -(\omega_g + \beta y)\mathbf{u}_g - f_0 \mathbf{u}_a - \partial_y \mathbf{k}_g, \quad (9b) $$

$$ \partial_t \mathbf{h} = -\partial_x(\mathbf{u}_g \mathbf{h}) - \partial_y(\mathbf{v}_g \mathbf{h}) - \mathbf{H}(\partial_x \mathbf{u}_a + \partial_y \mathbf{v}_a). \quad (9c) $$

In this form, the QG equations still appear similar to the RSW equations, but this is misleading as there is only one degree of freedom per layer. The two components of the geostrophic velocity are constrained by **h** via (7a) and (3a). Another difficulty with this system is that, although well-posed, it is highly implicit. There is no straightforward method to determine the ageostrophic velocity. In practice, and particularly for numerical models, the equations are rewritten in a more explicit form. This form reveals a hidden variable of the system that is the potential vorticity (PV):

$$ \mathbf{q} = \beta y + \omega_g - f_0 \frac{\mathbf{h}}{\mathbf{H}}. \quad (10) $$

## 2.3 Multi-layer QG equations

The multi-layer QG equations are obtained upon by differentiating (10) with respect to time and utilizing (9). The resulting equations are

$$ \partial_t \mathbf{q} = -\partial_x(\mathbf{u}_g \mathbf{q}) - \partial_y(\mathbf{v}_g \mathbf{q}), \quad (11a) $$

$$ \Delta \mathbf{p} - f_0^2 A \mathbf{p} = f_0 \mathbf{q} - f_0 \beta y, \quad (11b) $$

$$ -f_0 \mathbf{u}_g = \partial_y \mathbf{p}, \quad f_0 \mathbf{v}_g = \partial_x \mathbf{p}, \quad (11c) $$

where $\Delta = \partial_{xx}^2 + \partial_{yy}^2$ denotes the horizontal Laplacian operator, and $A$ represents the vertical discretization of the stretching operator introduced in Equation (3d). The PV $\mathbf{q}$ is the sole prognostic variable in (11a). All other model variables are derived from $\mathbf{q}$ via the diagnostic relations (11b) for $\mathbf{p}$ and (11c) for $\mathbf{u}$ and $\mathbf{v}$. The QG model is thus expressed in terms of the variables $(\mathbf{q}, \mathbf{p})$, with $\mathbf{p}/f_0$ being the vector of streamfunction for each layer.

Numerical implementations of the QG model typically rely on this system of equations. However, this system differs significantly from the RSW system of equations. One notable aspect is the absence of the ageostrophic velocities from the RSW system of equations. This absence does not imply that $\mathbf{u}_a$ and $\mathbf{v}_a$ vanish; rather, it means they can be disregarded if the focus is solely on the evolution of $\mathbf{p}$ and $\mathbf{q}$. The substantial difference between the QG and RSW model equations complicates direct comparisons between them. Consequently, numerical versions of the QG and RSW models generally have distinct implementations, with differing numerical methods in all aspects. In the following sections, we demonstrate how to reestablish the similarity between the two models.

# 3 Quasi-geostrophic model as projected rotating shallow-water

In this section, we propose a straightforward formulation of the QG equation using the horizontal velocity and layer thickness $(\mathbf{u}, \mathbf{v}, \mathbf{h})$ as state variables, along with the non-linear RSW operator $F$ defined in (6b), and a linear projection operator $P$. To our knowledge, this projection relation has been largely overlooked, with one notable exception being Charve (2004).

## 3.1 Projection formulation

To derive this projection, we start from the QG equations (8), and isolate the part controlled by $F$:

$$ \partial_t \begin{pmatrix} \mathbf{u}_g \\ \mathbf{v}_g \\ \mathbf{h} \end{pmatrix} = F \begin{pmatrix} \mathbf{u}_g \\ \mathbf{v}_g \\ \mathbf{h} \end{pmatrix} + \begin{pmatrix} f_0 \mathbf{v}_a \\ -f_0 \mathbf{u}_a \\ -\mathbf{H}(\partial_x \mathbf{u}_a + \partial_y \mathbf{v}_a) \end{pmatrix}, \quad (12) $$
211

212

213

214

where several terms cancel due to the geostrophic balance. In this form, the QG model can be interpreted as a RSW model applied to the geostrophic component and forced by an ageostrophic source term, represented by the second term of the right-hand side.

Let us now introduce the PV linear operator $Q$:

215

$$Q \begin{pmatrix} \mathbf{u} \\ \mathbf{v} \\ \mathbf{h} \end{pmatrix} = \partial_x \mathbf{v} - \partial_y \mathbf{u} - f_0 \frac{\mathbf{h}}{\mathbf{H}}. \qquad (13a)$$

216

This operator relates to the PV by

217

$$Q \begin{pmatrix} \mathbf{u}_g \\ \mathbf{v}_g \\ \mathbf{h} \end{pmatrix} = \mathbf{q} - \beta y, \qquad (13b)$$

218

and ensures that the contribution of the ageostrophic source term in (12) vanishes:

219

$$Q \begin{pmatrix} f_0 \mathbf{v}_a \\ -f_0 \mathbf{u}_a \\ -\mathbf{H}(\partial_x \mathbf{u}_a + \partial_y \mathbf{v}_a) \end{pmatrix} = 0. \qquad (13c)$$

220

Since $Q$ is independent of time, it commutes with the time derivative. Therefore,

221

$$\partial_t Q \begin{pmatrix} \mathbf{u}_g \\ \mathbf{v}_g \\ \mathbf{h} \end{pmatrix} = Q \circ F \begin{pmatrix} \mathbf{u}_g \\ \mathbf{v}_g \\ \mathbf{h} \end{pmatrix}. \qquad (14)$$

222

223

224

This equation is a reformulation of (11a), representing the conservation of PV. As $Q$ reduces the state vector from three to one variable, we need a reverse operation to retrieve the three variables $(\mathbf{u}_g, \mathbf{v}_g, \mathbf{h})$. This is accomplished by introducing the geostrophic operator $G$:

225

$$G(\mathbf{p}) = \begin{pmatrix} -\partial_y \mathbf{p} \\ \partial_x \mathbf{p} \\ f_0 \mathbf{H} A \mathbf{p} \end{pmatrix}. \qquad (15a)$$

226

227

Thus, $G$ is a linear operator with uniform and time-independent coefficients. The composition of the two operators $Q$ and $G$ yields:

228

$$(Q \circ G)(\mathbf{p}) = \Delta \mathbf{p} - f_0^2 A \mathbf{p} \qquad (15b)$$

229

230

which is the QG elliptic operator (11b) relating $\mathbf{p}$ to $\mathbf{q}$. This operator is invertible (Dutton, 1974; Bourgeois & Beale, 1994), hence $(Q \circ G)^{-1}$ is well-defined.

231

Finally, we introduce the QG projection operator $P$:

232

$$P = G \circ (Q \circ G)^{-1} \circ Q. \qquad (16)$$

233

$P$ is a projection because $P \circ P = P$, as demonstrated by

234

$$P \circ P = G \circ (Q \circ G)^{-1} \circ (Q \circ G) \circ (Q \circ G)^{-1} \circ Q = P. \qquad (17)$$

235

236

By construction, $P$ preserves the geostrophic state $\mathbf{X}_g = (\mathbf{u}_g, \mathbf{v}_g, \mathbf{h})^T$, i.e., $P(\mathbf{X}_g) = \mathbf{X}_g$. Indeed, since $\mathbf{X}_g$ is geostrophic, we have

237

$$Q(\mathbf{X}_g) = \mathbf{q} - \beta y, \qquad (18)$$

238

$$(Q \circ G)^{-1}(\mathbf{q} - \beta y) = \mathbf{p}/f_0, \qquad (19)$$

239

$$G(\mathbf{p}/f_0) = \mathbf{X}_g. \qquad (20)$$

240

241

Applying $P$ to Equation (12), and noting that $P$ commutes with $\partial_t$, we can formulate the multi-layer QG model as

242

$$\partial_t \mathbf{X}_g = P \circ F(\mathbf{X}_g). \qquad (21)$$
243

244

245

246

This form differs from the RSW model (6a) by the additional projection operator $P$ acting on $F$. By projecting the RSW tendency onto the geostrophic manifold, $P$ ensures that the QG state remains in geostrophic balance. Alternatively, the QG model can be viewed as evolving under the action of the RSW operator with the ageostrophic tendency removed.

247

248

249

250

251

252

253

254

To fully implement the projection method, it is necessary to specify the lateral boundary condition for Equation (15b). The appropriate condition, consistent with the QG model, is the Dirichlet boundary condition, where the pressure remains constant along each connected boundary. Given that the projection is applied to the model tendency rather than the model state, and that this work focuses on simply connected domains, the constant term is zero. In cases involving multiply connected domains, the constant can be explicitly evaluated as proposed by J. C. McWilliams (1977); however, this is beyond the scope of the present paper.

255

## 3.2 Diagnostic variables

256

257

With this formulation, the ageostrophic velocity ($u_a, v_a$) can be expressed in a simple closed form. Using equations (12) and (21), we have

258

$$ \begin{pmatrix} f_0 v_a \\ -f_0 u_a \\ -\mathbf{H}(\partial_x u_a + \partial_y v_a) \end{pmatrix} = P \circ F \begin{pmatrix} \mathbf{u}_g \\ \mathbf{v}_g \\ \mathbf{h} \end{pmatrix} - F \begin{pmatrix} \mathbf{u}_g \\ \mathbf{v}_g \\ \mathbf{h} \end{pmatrix}. \quad (22) $$

259

260

In a numerical model, this requires only a basic difference between the QG and RSW tendencies, with no costly computations.

261

262

263

264

265

266

267

268

269

Using this projection formulation, the QG system can be expressed with the same prognostic variables as the RSW system, namely the horizontal velocity ($u, v$) and the layer thickness $h$. This reestablishes the proximity between these two sets of equations in the ocean model hierarchy, leading to four practical implications. First, given a RSW discretization, one can implement the projection to derive the companion QG model for this specific discretization. Second, this formulation provides access to the ageostrophic velocity hidden in the QG equations at the cost of a simple subtraction. This enables an a posteriori diagnostic of the validity of the QG scaling. For a given geostrophic state ($u_g, v_g, h$), the corresponding ageostrophic velocity ($u_a, v_a$) can be computed and checked against:

270

$$ \frac{\mathbf{k}_a}{\mathbf{k}_g} \sim \mathrm{Ro}^2, \quad \frac{\mathbf{h}}{\mathbf{H}} \sim \mathrm{Ro}, \quad (23) $$

271

272

where $\mathbf{k}_a = (|\mathbf{u}_a|^2 + |\mathbf{v}_a|^2)/2$ and $\mathbf{k}_g = (|\mathbf{u}_g|^2 + |\mathbf{v}_g|^2)/2$ are the ageostrophic and the geostrophic kinetic energies, respectively.

273

274

275

276

277

278

Third, using the same variables ($u, v, h$) facilitates the development of eddy parameterizations, such as those in Bachman (2019); Li et al. (2023), which work similarly for both the RSW and the QG models. Finally, this formulation may offer new insights and approaches for the mathematical analysis of the QG equations, particularly regarding the extension of recent results on the well-posedness properties for a stochastic RSW model (Crisan & Lang, 2023) by examining the effect of the QG projector on these properties.

279

280

281

282

283

284

285

286

This formulation also presents an interesting analogy with the Leray projector formulation of the Navier-Stokes equations, where the Leray projector enforces the incompressible constraint by filtering out compressible sound waves through an elliptic Poisson equation for the pressure. Similarly, the proposed QG projector enforces the QG balance by filtering out fast gravity waves through an elliptic Helmholtz equation for the PV. It is noteworthy that the QG balance is not a constraint added to the RSW model; the Hamiltonian structure of the QG model (Holm & Zeitlin, 1998) is distinct from that of the RSW model. Specifically, contrary to the pressure, the PV is not a Lagrange multiplier enforcing a constraint.
# 4 Numerical experiments

This section presents the numerical experiments carried out to validate and explore the behavior of our new formulation. We first describe the implementation of the RSW solver and QG projector, then detail the specific experiments designed to test vortex shear instability, the idealized double-gyre circulation, and a simplified North Atlantic configuration. The first two experiments are analyzed to compare the solutions obtained from the QG and RSW models under varying conditions and parameters.

## 4.1 Numerical implementations

To implement our formulation, we require a RSW solver and a QG projector. We utilize the RSW solver developed by Roullet and Gaillard (2022), which is available at https://github.com/pvthinker/pyRSW. We re-implemented this solver in PyTorch to enable seamless GPU acceleration, verifying that our implementation reproduces their original results up to numerical precision. The key element of the discretization is a fifth-order Weighted Essentially Non-Oscillatory (WENO) upwinded reconstruction for both the mass flux and the nonlinear vortex-force term. This approach provides sufficient numerical dissipation, eliminating the need for an ad-hoc hyperviscous diffusion, while ensuring good material conservation of PV (Roullet & Gaillard, 2022).

The QG projector requires an elliptic solver. Thiry et al. (2024) released an efficient Python-PyTorch implementation of the QG equations. They solve the QG elliptic equation using discrete sine transforms implemented with PyTorch's Fast Fourier Transform (FFT), which leverages highly optimized MKL FFT on Intel CPUs and cuFFT on Nvidia GPUs. We employ their elliptic solver to solve the QG elliptic equation (11b). Once the projector is available, the QG solver is implemented by simply adding the projection step to the RSW solver. Numerically, the two solvers differ by only one line of code – the projection step. All other elements, including variable staggering, time discretization, and core RSW equations, are identical.

To summarize, the QG model is advanced in time as follows: (1) compute the RSW tendencies, (2) apply the projection operator $P$ to obtain the QG tendencies, and (3) update the state variables. For the RSW model, step (2) is omitted, and only steps (1) and (3) are performed.

Consequently, we achieve a PyTorch implementation (see Supporting Information files) that is concise (approximately 800 lines of code), remains true to the equations, and implements both the multi-layer QG and RSW equations using the same state variables ($u, v, h$).

## 4.2 Vortex shear instability

To validate our formulation, we first study a vortex shear instability and compare its evolution in both the QG and the RSW models as $Ro$ increases for $Bu = 1$. The initial state is a perfectly shielded vortex with a core of uniform vorticity $\omega_1$ surrounded by a ring of opposite sign uniform vorticity $\omega_2$. The system involves two lengths: the core radius $r_0$ and the vortex outer radius $r_1$. The ratio $\omega_2/\omega_1$ is such that the total circulation vanishes. This system is unstable and leads to the formation of multipoles (Morel & Carton, 1994). The number of poles depends on the ratio of the vortex radius to the core radius. We focus here on a tripole formation case with $r_1/r_0 = 1.4$. To promote the growth of the most unstable mode, we add a small mode azimuthal perturbation.

The experiments are set up in dimensional form with a square domain of size $L_x \times L_y = 100$ km $\times$ 100 km on a f-plane. There is a single layer of fluid with thickness $H = 1$ km. We assume no-flow and free-slip boundary conditions. The acceleration of gravity is set to $g = 10 \text{ m s}^{-2}$. The vortex core has a radius $r_0 = 10$ km and a positive vorticity; the vortex radius is $r_1 = 14$ km. We run the simulations on a 512 $\times$ 512 grid, providing a 200m spatial
resolution. The Coriolis parameter is chosen such that $Bu = 1$, i.e., $f_0 = \sqrt{gH}/r_0$. The vorticity $\omega_1$ is set indirectly via $u_{max}$, the maximum velocity of the initial condition. We define the Rossby number as $Ro = u_{max}/(f_0r_0)$. We impose $Ro$ and deduce $u_{max}$. We have tested four cases: $Ro \in \{0.01, 0.05, 0.1, 0.5\}$. In all cases, we apply the QG projector to the initial state before starting the time integration. Consequently, the initial states are in geostrophic balance and differ only by a scaling factor. As $Ro$ increases, the instability develops faster in dimensional time. Therefore, we present the results in the rescaled eddy-turnover time $\tau$, defined as the inverse of the $l^2$-norm of the initial vorticity, i.e., $\tau = 1/\|\omega\|$. We integrate the simulations for $T = 10\,\tau$.

Figure 2 shows the evolution of the initial state for the two extreme cases, $Ro = 0.01$ and $Ro = 0.5$. In the $Ro = 0.01$ case, the differences between the QG and the RSW solutions are not perceptible to the naked eye. This is consistent with the QG scaling and the fact that the QG model is an asymptotic limit of the RSW model. In the $Ro = 0.5$ case, the differences are of order one. In the RSW model, the vortex spins faster, and the filaments of negative vorticity are stabilized. The smoothness of the RSW solution at $t = 10\,\tau$ might be surprising as we expect some gravity waves generated during the fast cyclo-geostrophic adjustment of the initial state. These gravity waves are present but at $t = 10\,\tau$, they have bounced back and forth several times along the boundary and are strongly scattered. Interestingly, the QG solution is exactly the same as in the $Ro = 0.01$ case. This is expected since the QG equations are scaling invariant, i.e., the evolution is invariant under multiplication by a constant. However, it is also quite remarkable to recover this property in the numerical solutions because the QG solver relies on the full RSW right-hand side. Another symmetry of the QG model is parity invariance. The solution should be invariant under a sign change of the vorticity. In other words, cyclones and anticyclones follow the same evolution, up to a change in rotation. This is not at all the case for the RSW model. By flipping the sign of the initial vorticity, we verified that the solutions satisfy this property (not shown). This means that the QG projector behaves well as it restores two invariances, scaling and parity, that are absent in the RSW solver.

To assess the influence of $Ro$ on the time evolution, we define the normalized difference $\delta = 2||\omega_{qg} - \omega_{sw}||/(||\omega_{qg}|| + ||\omega_{sw}||)$. Figure 3 shows $\delta(t)$ for the four $Ro$ cases. The oscillations during the $[0, 4\,\tau]$ period are due to the gravity waves in the RSW case resulting from the imperfect initial balance. Interestingly, the shortening of the oscillation period as $Ro$ decreases is due to the time rescaling by $\tau$. In dimensional time, the periods are the same. A practical consequence is that the RSW experiment requires more time steps to reach $t = 10\,\tau$ as $Ro$ decreases, whereas for the QG experiment, the number of time steps is constant. After $t = 4\,\tau$, the differences are dominated by the shear instability developing on the vortex. The small difference seen in the snapshots in the $Ro = 0.01$ case has actually reached a plateau, meaning the QG solution closely follows the RSW one. This confirms that the QG model is a good simplified model when the scaling assumption holds. On a longer time scale and with chaotic vortex dynamics, the solutions would diverge, but in such a simple setup, the solutions remain close. The time evolution of $\delta$ is similar for all $Ro$ except $Ro = 0.5$, where $\delta$ saturates at one, the maximum value of this metric. In this case, this metric suggests that the two models predict a completely different solution. Looking back at the snapshots, this seems exaggerated. The order one difference is mostly due to the difference in timing, not the difference in pattern. Depending on the purpose, the QG solution might still be of interest as it still captures the main phenomenology. Note that the QG solution could be made closer to the RSW solution if the two models were started with different initial states: a cyclo-geostrophic balance for the RSW model and the associated projected state for the QG model. The initial state would thus depend on $Ro$. For this illustrative experiment, we preferred to stick to the same initial state, up to a scaling constant.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Vortex shear instability solved using QG and SW for $Bu = 1$. (Top to bottom) Initial relative vorticity for $Ro = 0.01$ and $Ro = 0.5$. Final relative vorticities and their differences for QG and SW models with $Ro = 0.01$. Final relative vorticities and their differences for QG and SW models with $Ro = 0.5$. Units are in $s^{-1}$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Time evolution of the normalized differences of relative vorticity between the QG model and the RSW model with $Bu = 1$ and different $Ro$.</figcaption></figure>

<table><caption>**Table 1.** Parameters of the idealized double-gyre configuration. There no viscosity nor diffusion coefficient as the grid-scale dissipation is implicitly handled by the upwinded WENO reconstruction of the mass flux and the vortex-force (see Roullet & Gaillard, 2022).</caption><thead><tr><th>Parameter</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>$L_x, L_y$</td><td>5120, 5120 km</td><td>Domain size</td></tr><tr><td>$H_k$</td><td>400, 1100, 2600 m</td><td>Mean layer thickness</td></tr><tr><td>$g'_k$</td><td>9.81, 0.025, 0.0125 m s<sup>-2</sup></td><td>Reduced gravity</td></tr><tr><td>$\gamma$</td><td>3.6 10<sup>-8</sup> s<sup>-1</sup></td><td>Bottom drag coefficient</td></tr><tr><td>$\tau_0$</td><td>0.08 N m<sup>-2</sup></td><td>Wind stress magnitude</td></tr><tr><td>$\rho$</td><td>1000 kg m<sup>-3</sup></td><td>Ocean density</td></tr><tr><td>$f_0$</td><td>9.375 10<sup>-5</sup> s<sup>-1</sup></td><td>Mean Coriolis parameter</td></tr><tr><td>$\beta$</td><td>1.754 10<sup>-11</sup> (m s)<sup>-1</sup></td><td>Coriolis parameter gradient</td></tr><tr><td>$L_d$</td><td>41, 25 km</td><td>Baroclinic Rossby radii</td></tr><tr><td>$n_x, n_y$</td><td>256, 256</td><td>Grid size</td></tr><tr><td>$dt$</td><td>4000 s</td><td>Time-step</td></tr></tbody></table>

## 4.3 Double-gyre circulation

To explore a richer phenomenology, we have tested our new formulation on a classical oceanic test case, the idealized wind-forced double-gyre circulation. The domain is a non-periodic square ocean basin with $N = 3$ layers. We assume free-slip boundary conditions on each boundary. A stationary and symmetric wind stress ($\tau_x, \tau_y$) is applied in the top layer, with $\tau_x = -(\tau_0/\rho_0) \cos(2\pi y/L_y)$ and $\tau_y = 0$. Additionally, a linear drag with drag coefficient $\gamma$ is applied in the bottom layer. The parameter values are given in Table 1.

We study this configuration in an eddy-permitting resolution of 20 km, meaning that the spatial resolution (20 km) is half of the largest baroclinic Rossby radius (41 km). We expect QG and RSW simulations to produce strong western boundary currents converging to the middle of the western boundary and an eastward jet departing from the middle of the western boundary. Starting from the rest state, this configuration requires about 100 years to spin up and achieve converged statistics (Hogg et al., 2005; Simonnet, 2005).

In Figure 4, we present snapshots of key quantities of the upper layer after 100 years of integration. The solution exhibits expected properties of this setup: a double-gyre circulation (Fig. 4e) separated by a meandering eastward jet (Fig. 4a,d) emanating from the western boundary, two strong and narrow western boundary currents (Fig. 4b) feeding the eastward
jet, mesoscale turbulence throughout, especially near the jet (Fig. 4c,d), and Rossby waves propagation (Fig. 4e,f). The local Rossby number, defined as $\omega/f_0$ peaks at 0.2, which is significant but not large enough to dismiss the solution. Interestingly, the layer thickness perturbation $h_1$ is substantial compared to the reference depth $H_1 = 400$ m.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** QG upper-layer snapshot after 100 years of spin-up from the rest state. KE, geos., and rel. stand for kinetic energy, geostrophic, and relative, respectively.</figcaption></figure>

To complement this snapshot, we present statistics estimated over the years (100-200) after the initial state, using one snapshot every 10 days. We decompose the geostrophic kinetic energy into its mean and eddy components, and compute the ratio $k_a/k_g$. The results are shown in Figure 5 for the upper layer. The statistics confirm the presence of the strong western boundary currents (Fig. 5a) and a strong eastward jet reaching the middle of the domain (Fig. 5a,d). A remarkable feature is the symmetry of all quantities with respect to the central latitude, due to the symmetry of the forcing and the cyclone/anticyclone symmetry of the QG model. This symmetry prevents the model from breaking the symmetry of the forcing. From the ratio $k_a/k_g$ and $h/H$ (Fig. 5c,f), we test the validity of the QG scaling assumptions (Eq. 23). The color scale is adjusted so that the white intermediate color corresponds to $Ro = 0.1$. Red areas indicate where the QG scaling assumptions are not respected. Counter-intuitively, the central jet region, where kinetic energy is the highest, shows the best QG scaling. The worst regions are the gyre centers due to large thickness deviations and the boundaries due to large ageostrophic velocities.

Finally, we compare the QG solution with its companion RSW one. Due to the free surface and the presence of fast barotropic gravity waves, the RSW solver requires a much smaller time step, typically 200 times smaller. This factor corresponds to the ratio $\sqrt{g\bar{H}}/\max(u)$. Integrating from the rest state over 200 years would be computationally intensive. A compromise would be to use a barotropic-baroclinic time splitting (Higdon & De Szoeke, 1997) or implicit stepping of the free surface (Roullet & Madec, 2000), but this is beyond the scope of the present study. Instead, we ran the RSW solver starting from year 200 of the QG solution and integrate it over 2 years. Figure 6 compares QG and RSW on a snapshot of vorticity in the upper layer. The central jet has a southward component, the two gyres are no longer symmetric, and mesoscale turbulence has intensified in the Northern gyre and weakened in the Southern gyre. Applying the QG projector on the RSW state results in a state that remains very close (Fig. 6c) yet with damped fluctuations. The QG
<figure><img src="image_7.png" /><figcaption>**Figure 5.** QG upper-layer statistics over 100 years after 100 years spin-up from rest state.</figcaption></figure>

projector tends to dampen the short scales, but since it is applied on the RSW model tendency, it does not affect the QG state, as seen from the comparison between Fig. 6a and Fig. 6b. The closeness suggests that while the QG solution might be locally tangent to the RSW solution, long-term integration results in noticeable differences.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (Left) QG upper-layer relative vorticity after 100 years of spin-up, (middle) RSW upper-layer relative vorticity, and (right) projected RSW vorticity after 2 additional years of spin-up. Units are in $f_0$.</figcaption></figure>

## 4.4 Simplified North Atlantic configuration

To demonstrate the applicability of the proposed projection formulation to more complex ocean configurations, we extend the previous idealized double-gyre circulation to the North Atlantic basin, located at 9°N–48°N and 98°W–4°W. The ocean boundaries are set at a depth of 250 meters to exclude the continental shelf from the simulation. In this study, we simplify the realistic North Atlantic configuration by not considering inflows and outflows across the northern and southern open boundaries (Marchesiello et al., 2001; Blayo & Debreu, 2005), which will be investigated in subsequent studies. Here, we assume a flat bottom to be consistent with the formulation presented in the previous sections. However, the presented formulation could be extended to include bottom topography in the future. The
ocean surface boundary condition is imposed using the Hellerman and Rosenstein (1983) monthly wind stress climatology, which has been widely used in several other reference simulations of the Gulf Stream (Blayo et al., 1994; Hurlburt & Hogan, 2000).

The simulation runs on a horizontal grid of 1024 × 512 points, corresponding to a spatial resolution of 8.5 km (~1/11°), with 3 vertical layers using the same reference thickness and reduced gravity as specified in Table 1. The time step is set to 2000 s, and a partial free-slip condition is applied on the ocean boundaries with the coefficient set to 0.6 (0 for no-slip and 1 for free-slip). To solve the elliptic equation on non-rectangular geometries for the QG projector, we extend the fast discrete sine transform spectral solver using the capacitance matrix method (Thiry et al. (2024), see also Blayo and LeProvost (1992) for the original technique).

Figure 7 illustrates the instantaneous prognostic and diagnostic variables of the projected QG model after 40 years of spin-up simulation. The surface relative vorticity in Fig. 7a reveals the meandering and energetic western boundary current, along with the turbulent eddying structures throughout most of the basin. The sea surface height in Fig. 7b illustrates the gyre circulation structure in this configuration and could potentially be assimilated with available observations. Fig. 7c and Fig. 7d demonstrate the surface current speed ($\sqrt{u_1^2 + v_1^2}$) for the prognostic geostrophic motions and the diagnostic ageostrophic motions, respectively. The latter is at least five times smaller than the former across the entire basin. This instantaneous snapshot notably shows that the ageostrophic velocity component is extremely small in the region going from 25°N to 30°N, which approximately corresponds to the range of latitudes for the state of Florida. However, the ageostrophic motions remain significant along the Gulf Stream current.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** QG upper-layer snapshot after 40 years of spin-up from the rest state.</figcaption></figure>

We also investigate the temporal statistics of the projected QG model over 20 years following a 40-year spin-up. For instance, Figure 8 illustrates the surface MKE and EKE densities for both geostrophic and ageostrophic motions. Figures 8a and 8b reveal that EKE predominates over MKE across the basin, particularly along the Gulf Stream, highlighting the significant influence of mesoscale eddies on flow variability. Figures 8c and 8d show the ratio of ageostrophic to geostrophic MKE and EKE, confirming the previous instantaneous
finding: geostrophic balance is markedly dominant in the North Atlantic Subtropical Gyre, validating the QG scaling in that region.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** QG upper-layer statistics over 20 years after 40 years spin-up from rest state.</figcaption></figure>

# 5 Conclusions

In this paper, we have demonstrated that the QG model can be formulated as a projected RSW model by applying the QG projector $P$ (defined in Equation 16) to the RSW tendency. This formulation allows the QG model to utilize the same variables $(u, v, h)$ as the RSW model, rather than the more conventional $(p, q)$ variables. This unified approach enables us to enhance the similarities and differences between the two models and facilitates a coherent integration of both models within the same numerical framework. In contrast to earlier approaches based on similar concepts (Leith, 1980; Salmon, 1998; Saujani & Shepherd, 2006), the projection formulation we propose here does not depend on the shape of the boundary. As demonstrated, it can be easily implemented for basins with general coastal geometries.

We have validated this approach through a vortex shear instability test, demonstrating that the resulting QG model retains essential symmetry properties, such as scaling and parity, which are not present in the RSW model. Furthermore, we have shown that a QG model implemented in this manner reproduces the expected characteristics of a double-gyre experiment. An immediate benefit of this approach is the ability to directly compare RSW and QG solutions. This straightforward formulation paves the way for further investigations into the differences between QG and RSW equations, such as studies on the stability of geostrophic equilibrium or spontaneous imbalance phenomena.

Additionally, it holds significant potential for data assimilation applications, enabling seamless switching between these nested RSW and QG models or even their combination. For instance, one could initially run a data assimilation algorithm, such as an ensemble Kalman filter (Evensen, 2003), using the projected QG model to capture geostrophic dynamics. Subsequently, the obtained solution could be refined by transitioning to the RSW model. This flexibility could enhance the accuracy and efficiency of data assimilation processes in oceanographic modeling and related fields.
504

## Open Research Section

505

506

The code files necessary to reproduce the results presented in this paper are publicly available in Thiry and Li (2024).

507

508

## Acknowledgments

509

510

511

512

513

514

515

The authors gratefully acknowledge the two anonymous reviewers for their insightful comments and suggestions, which helped to significantly improve the manuscript. We also wish to acknowledge the support of the ERC EU project 856408-STUOD, which enabled the research presented in this paper. Additionally, the authors extend their thanks to Anne-Laure Dalibard for her invaluable contribution in directing them to Frederic Charve's Thesis. Furthermore, the authors express their gratitude to Frederic Charve for engaging discussions.

516

## References

517

518

Bachman, S. (2019). The GM+E closure: a framework for coupling backscatter with the Gent and McWilliams parameterization. *Ocean Modelling*, **136**, 85–106.

519

520

Bachman, S., Fox-Kemper, B., & Pearson, B. (2017). A scale-aware subgrid model for quasi-geostrophic turbulence. *Journal of Geophysical Research: Oceans*, **122**, 1529–1554.

521

522

523

Bauer, W., Chandramouli, P., Chapron, B., Li, L., & Mémin, E. (2020). Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: A stochastic approach. *Journal of Physical Oceanography*, **50**, 983–1003.

524

525

Berloff, P. (2005). Random-forcing model of the mesoscale oceanic eddies. *Journal of Fluid Mechanics*, **529**, 71–95.

526

527

Blayo, E., & Debreu, L. (2005). Revisiting open boundary conditions from the point of view of characteristic variables. *Ocean Modelling*, **9**(3), 231–252.

528

529

530

Blayo, E., & LeProvost, C. (1992, 10). Performance of the capacitance matrix method for solving Helmholtz-type equations in ocean modelling. *Journal of Computational Physics*, **102**.

531

532

533

Blayo, E., Verron, J., & Molines, J. M. (1994). Assimilation of TOPEX/POSEIDON altimeter data into a circulation model of the North Atlantic. *Journal of Geophysical Research: Oceans*, **99**(C12), 24691–24705.

534

535

536

Bourgeois, A., & Beale, J. (1994). Validity of the quasigeostrophic model for large-scale flow in the atmosphere and ocean. *SIAM Journal on Mathematical Analysis*, **25**(4), 1023–1068.

537

538

Charve, F. (2004). *Etude de phénomènes dispersifs en mécanique des fluides géophysiques* (Unpublished doctoral dissertation). Palaiseau, Ecole polytechnique.

539

540

Crisan, D., & Lang, O. (2023). Well-posedness properties for a stochastic rotating shallow water model. *Journal of Dynamics and Differential Equations*, 1–31.

541

542

Cushman-Roisin, B., & Beckers, J. (2011). *Introduction to geophysical fluid dynamics: physical and numerical aspects*. Academic Press.

543

544

Dutton, J. (1974). The nonlinear quasi-geostrophic equation: Existence and uniqueness of solutions on a bounded domain. *Journal of Atmospheric Sciences*, **31**(2), 422–433.

545

546

Evensen, G. (2003). The ensemble Kalman filter: Theoretical formulation and practical implementation. *Ocean Dynamics*, **53**, 343–367.

547

548

Grooms, I., Majda, A., & Smith, K. (2015). Stochastic superparameterization in a quasi-geostrophic model of the Antarctic Circumpolar Current. *Ocean Modelling*, **85**, 1–15.

549

550

Hellerman, S., & Rosenstein, M. (1983). Normal monthly wind stress over the world ocean with error estimates. *Journal of Physical Oceanography*, **13**, 1093–1104.

551

552

Higdon, R., & De Szoeke, R. (1997). Barotropic-baroclinic time splitting for ocean circulation modeling. *Journal of Computational Physics*, **135**(1), 30–53.
553

554

555

Hogg, A., Killworth, P., Blundell, J., & Dewar, W. (2005). Mechanisms of decadal variability of the wind-driven ocean circulation. *Journal of Physical Oceanography*, 35(4), 512–531.

556

557

Holm, D. D., & Zeitlin, V. (1998). Hamilton's principle for quasigeostrophic motion. *Physics of fluids*, 10(4), 800–806.

558

559

Holton, J. (1973). An introduction to dynamic meteorology. *American Journal of Physics*, 41(5), 752–754.

560

561

562

Hurlburt, H. E., & Hogan, P. J. (2000). Impact of 1/8° to 1/64° resolution on Gulf Stream model–data comparisons in basin-scale subtropical Atlantic ocean models. *Dynamics of Atmospheres and Oceans*, 32(3), 283–329.

563

564

Jansen, M., & Held, I. (2014). Parameterizing subgrid-scale eddy effects using energetically consistent backscatter. *Ocean Modelling*, 80, 36–48.

565

566

Leith, C. E. (1980). Nonlinear normal mode initialization and quasi-geostrophic theory. *Journal of the Atmospheric Sciences*, 37, 958–968.

567

568

569

Li, L., Deremble, B., Lahaye, N., & Mémin, E. (2023). Stochastic data-driven parameterization of unresolved eddy effects in a baroclinic quasi-geostrophic model. *Journal of Advances in Modeling Earth Systems*, 15(2), e2022MS003297.

570

571

Marchesiello, P., McWilliams, J. C., & Shchepetkin, A. (2001). Open boundary conditions for long-term integration of regional oceanic models. *Ocean Modelling*, 3(1), 1–20.

572

573

Marshall, D., Maddison, J., & Berloff, P. (2012). A framework for parameterizing eddy potential vorticity fluxes. *Journal of Physical Oceanography*, 42(4), 539–557.

574

575

McWilliams, J. (2006). *Fundamentals of geophysical fluid dynamics*. Cambridge University Press.

576

577

McWilliams, J. C. (1977). A note on a consistent quasigeostrophic model in a multiply connected domain. *Dynamics of Atmospheres and Oceans*, 1(5), 427–441.

578

579

Morel, Y., & Carton, X. (1994). Multipolar vortices in two-dimensional incompressible flows. *Journal of Fluid Mechanics*, 267, 23–51.

580

581

582

Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., ... others (2019). Pytorch: An imperative style, high-performance deep learning library. *Advances in Neural Information Processing Systems*, 32.

583

Pedlosky, J. (2013). *Geophysical fluid dynamics*. Springer Science & Business Media.

584

585

586

Roullet, G., & Gaillard, T. (2022). A fast monotone discretization of the rotating shallow water equations. *Journal of Advances in Modeling Earth Systems*, 14(2), e2021MS002663.

587

588

589

Roullet, G., & Madec, G. (2000). Salt conservation, free surface, and varying levels: a new formulation for ocean general circulation models. *Journal of Geophysical Research: Oceans*, 105(C10), 23927–23942.

590

591

592

Salmon, R. (1998). *Lectures on geophysical fluid dynamics*. Oxford Academic.

Saujani, S., & Shepherd, T. G. (2006). A unified theory of balance in the extratropics. *Journal of Fluid Mechanics*, 569, 447–464.

593

594

Simonnet, E. (2005). Quantization of the low-frequency variability of the double-gyre circulation. *Journal of Physical Oceanography*, 35(11), 2268–2290.

595

596

597

Thiry, L., & Li, L. (2024). *A unified formulation of quasi-geostrophic and shallow water equations via projection*. Zenodo. Retrieved from https://doi.org/10.5281/zenodo.13372736 ([Software])

598

599

600

Thiry, L., Li, L., Roullet, G., & Mémin, E. (2024). MQGeometry-1.0: a multi-layer quasi-geostrophic solver on non-rectangular geometries. *Geoscientific Model Development*, 17(4), 1749–1764.

601

602

Uchida, T., Dersemble, B., & Popinet, S. (2022). Deterministic model of the eddy dynamics for a midlatitude ocean model. *Journal of Physical Oceanography*, 52(6), 1133–1154.

603

604

605

Vallis, G. (2017). *Atmospheric and oceanic fluid dynamics*. Cambridge University Press.

Zanna, L., Mana, P., Anstey, J., David, T., & Bolton, T. (2017). Scale-aware deterministic and stochastic parametrizations of eddy-mean flow interaction. *Ocean Modelling*, 111, 66–80.
Zeitlin, V. (2018). *Geophysical fluid dynamics: understanding (almost) everything with rotating shallow water models*. Oxford University Press.