# Stochastic data-driven parameterization of unresolved mesoscale eddies

Long Li, Bruno Deremble, Noé Lahaye, Etienne Mémin

► To cite this version:

Long Li, Bruno Deremble, Noé Lahaye, Etienne Mémin. Stochastic data-driven parameterization of unresolved mesoscale eddies. 2022. hal-03727820

HAL Id: hal-03727820

https://hal.science/hal-03727820v1

Preprint submitted on 19 Jul 2022

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Stochastic data-driven parameterization of unresolved mesoscale eddies

Long Li*

Bruno Deremble†

Noé Lahaye*

Etienne Mémin*

## Abstract

In this work, a stochastic representation based on a physical transport principle is proposed to account for mesoscale eddy effects on the large-scale oceanic circulation. This stochastic framework arises from a decomposition of the Lagrangian velocity into a smooth-in-time component and a highly oscillating noise term. One important characteristic of this random model is that it conserves the total energy of the resolved flow for any realization. Such an energy-preserving representation is successfully implemented in a well established multi-layered quasi-geostrophic dynamical core. The empirical spatial correlation of the unresolved noise is calibrated from the eddy-resolving simulation data. In particular, a stationary correction drift can be introduced in the noise through Girsanov transformation. This non intuitive term appears to be important in reproducing on a coarse mesh the eastwards jet of the wind-driven double-gyre circulation. In addition, a projection method has been proposed to constrain the noise living along the iso-surfaces of the vertical stratification. The resulting noise enables us to improve the intrinsic low-frequency variability of the large-scale current.

## 1 Introduction

Ocean mesoscale eddies contain a large proportion of energy and have an important impact on large-scale circulations. They are found everywhere and are particularly intensive in the western boundary currents such as the Gulf Stream and the Kuroshio. However, these eddies are generally not resolved in ocean general circulation models, since the deformation radius in the ocean is at most of the order of 100 km. In particular, the effects of the mesoscale eddies need to be properly parameterized in coarse-resolution ocean models.

The most successful parameterization of mesoscale effect is based on the eddy-induced advection scheme (Gent and McWilliams, 1990; Gent et al., 1995; Griffies, 1998), which consists in mimicking the impact of the baroclinic instability by flattening the isopycnal surfaces to transfer the available potential energy of the resolved flow to the eddy kinetic energy at the sub-grid scales. However, such scheme does not account for the backscattering of kinetic energy from small to large scales (Bachman, 2019). Classical eddy viscosity models are introduced in coarse models to mimic the mixing action of the unresolved scales. The associated energy dissipation mechanism is often represented by some functional operators (Leith, 1971; Griffies and Hallberg, 2000; Bachman et al., 2017) that depend on the resolved flow. These approaches seem to improve the effective resolution of eddy-permitting models. Nevertheless, they still lead to excessive decreases of the resolved kinetic energy in large-scale ocean models. The lack of variability due to energy missing is particularly detrimental in ensemble forecasting and data assimilation applications (Karspeck et al., 2013; Franzke et al., 2015).

<footnote>*Centre Inria Rennes – Bretagne Atlantique, France</footnote>

<footnote>†Institut des Géosciences de l'Environnement, Grenoble, France</footnote>
To address these shortcomings, stochastic parameterizations are becoming more and more popular and several schemes have been devised and studied for ocean circulation models. These schemes introduce energy backscattering across scales and enable to increase the internal ocean variability. For instance, Berloff (2005) developed a random-forcing model based on the dynamical decomposition of the flow into a large-scale component and an eddy components; Jansen and Held (2014) proposed a stochastic forcing in the form of a negative Laplacian viscosity to inject dissipated energy back to the resolved flow; Grooms and Majda (2014); Grooms et al. (2015) proposed a stochastic superparameterization scheme that includes stochastic Reynolds stress terms to backscatter kinetic energy; Porta Mana and Zanna (2014); Zanna et al. (2017) proposed an eddy parameterization based on a non-Newtonian stress which depends on the partially resolved scales and their variability.

We propose here a specific stochastic representation of the unresolved flow dynamics based on the location uncertainty (LU) framework (Mémin, 2014). This representation introduces a random unresolved flow component and relies on a stochastic transport principle. Similar approaches based on the same decomposition have also been recently proposed by Holm (2015); Cotter et al. (2019, 2020); Gugole and Franzke (2019). This kind of random models enable to consider systems that are less dissipative than their deterministic counterparts. Nevertheless, the ability of such models to represent faithfully the uncertainties associated to the actual unresolved small scales highly depends on the definition of the random component and on its evolution in time. Unsurprisingly, stationarity/time-varying and homogeneity/inhomogeneity characteristics of the unresolved flow component have strong influences on the numerical results. For instance, Bauer et al. (2020a) illustrates that the noise inhomogeneity induces a structuration of the large-scale flow reminiscent to the action of the vortex force associated to surface wave-induced Stokes drift; Bauer et al. (2020b) shows that the introduction of inhomogeneous noise into the barotropic quasi-geostrophic (QG) model enables us to reproduce accurately, on a coarse mesh, the high order statistics of the eddy-resolving data; a stochastic shallow water model preserving the resolved total energy has been proposed by Brecht et al. (2021) and the results show that the stochastic parameterization provides a good trade-off between model error representation and ensemble spread.

In this work, we investigate the LU formulation for the baroclinic QG model (Hogg et al., 2004) within an idealized double-gyre circulation configuration. In particular, we focus on the reproduction of the eastwards jet as well as the prediction of the low-frequency variability for the proposed random model on a coarse resolution grid, for cases in which the baroclinic instability is not resolved. Another, important aspect, concerns the ability to include in the noise representation a stationary drift component associated to the temporal mean of the high-resolution fluctuations. As shown in this paper such stationary drift can be elegantly introduced in the noise through Girsanov transformation. In addition, a projection method has been proposed to update in time the noise along the iso-surfaces of the stratification.

The remainder of this paper is structured as follows. Section 2 describes the framework of LU and presents the derived stochastic QG model. Section 3 explains the discretization of the proposed QG model and the parameterization method for the unresolved noise. Section 4 discusses the numerical results with some statistical diagnosis and energetic analysis. In Section 5 we draw some conclusions and provide an outlook for future works. In the Appendices we demonstrate the energy conservation of the random system and explore the energy conversion between the ensemble-mean and the ensemble-eddy components.
# 2 Continuous models

In this section, we first review the general setting of Location Uncertainty (LU) (Mémin, 2014), we then present the resulting stochastic Boussinesq quasi-geostrophic (QG) model for a continuously stratified ocean, finally we complete the random formulation by including a time-correlated drift into the unresolved flow component through the Girsanov transformation.

## 2.1 Stochastic flow

The evolution of Lagrangian particle trajectory ($X_t$) under LU is described by the following stochastic differential equation (SDE):

$$dX_t(x) = v(X_t(x), t) dt + \sigma(X_t(x), t) dB_t, \quad X_0(x) = x \in \mathcal{D}, \quad (2.1)$$

where $v$ denotes the time-smooth resolved velocity that is both spatially and temporally correlated, $\sigma dB_t$ stands for the fast oscillating unresolved flow component (also called noise in the following) that is correlated in space yet uncorrelated in time, and $\mathcal{D} \subset \mathbb{R}^d$ ($d = 2$ or 3) is a bounded spatial domain.

Mathematically, {$B_t$}$_{0 \le t \le T}$ is an $I_d$-cylindrical Brownian motion (Da Prato and Zabczyk, 2014) on a filtered probability space $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{0 \le t \le T}, \mathbb{P})$<fnref n="1" /> and takes values in the Hilbert space $H := (L^2(\mathcal{D}))^d$. Informally, one can consider its time-derivative as a space-time “white noise” (in a distribution sense) with independent components ($B_t^i$)$_{i=1,...,d}$.

The spatial structure of the unresolved flow component is modeled by the correlation operator, $\sigma$, which is not necessarily deterministic. In the most general case it is random and time-dependent with some regularity conditions. More precisely, for each $(\omega, t) \in \Omega \times [0, T]$, $\sigma(\cdot, t)[\bullet]$ is assumed to be a Hilbert-Schmidt integral operator on $H$ with a bounded matrix kernel $\check{\sigma} = (\check{\sigma}_{ij})_{i,j=1,...,d}$ such that

$$\sigma(x, t) f = \int_{\mathcal{D}} \check{\sigma}(x, y, t) f(y) dy, \quad f \in H, \quad x \in \mathcal{D}. \quad (2.2)$$

The composition of $\sigma(\cdot, t)[\bullet]$ and its adjoint $\sigma^*(\cdot, t)[\bullet]$ is a trace class operator on $H$ and admits eigenfunctions $\xi_n(\cdot, t)$ with eigenvalues $\lambda_n(t)$ satisfying $\sum_{n \in \mathbb{N}} \lambda_n(t) < +\infty$. Then, the noise can be equally defined by the following spectral decomposition:

$$\sigma(x, t) dB_t = \sum_{n \in \mathbb{N}} \lambda_n^{1/2}(t) \xi_n(x, t) d\beta_n(t), \quad (2.3)$$

where $\beta_n$ are independent (one-dimensional) standard Brownian motions.

In addition, we assume that the operator-space-valued process {$\sigma(\cdot, t)[\bullet]$}$_{0 \le t \le T}$ is stochastically integrable, i.e. $\mathbb{P}[\int_0^T \sum_{n \in \mathbb{N}} \lambda_n(t) dt < +\infty] = 1$. As such, the stochastic integral {$\int_0^t \sigma(\cdot, s) dB_s$}$_{0 \le t \le T}$ is a $H$-valued Gaussian process of zero mean and of bounded (global) variance under the probability measure $\mathbb{P}$:

$$\mathbb{E}_{\mathbb{P}} \left[ \int_{0}^{t} \sigma(\cdot, s) dB_{s} \right] = 0, \quad \mathbb{E}_{\mathbb{P}} \left[ \left\| \int_{0}^{t} \sigma(\cdot, s) dB_{s} \right\|_{H}^{2} \right] < +\infty. \quad (2.4)$$

<footnote n="1">¹$\Omega$ is a set of samples, $\mathcal{F}$ is a $\sigma$-field (a collection of subsets of $\Omega$), {$\mathcal{F}_t$}$_{0 \le t \le T}$ is a filtration (a family of sub $\sigma$-fields of $\mathcal{F}$ that are continuously indexed in time and ordered non-decreasingly), and $\mathbb{P}$ is a probability measure.</footnote>
Moreover, the point-wise ($x \in \mathcal{D}$) and path-wise ($\omega \in \Omega$) strength of the unresolved flow component at each instant ($t \in [0, T]$) is measured by the matrix kernel of the composite operator $\sigma\sigma^*$, and denoted by $a(x, t)$, namely

$$a(x,t) := \int_{\mathcal{D}} \check{\sigma}(x, y, t) \check{\sigma}^T(y, x, t) dy = \sum_{n \in \mathbb{N}} \lambda_n(t) (\xi_n \xi_n^T)(x, t). \quad (2.5)$$

As shown in Bauer et al. (2020a), the process $\int_0^t a(x, s) ds$ actually corresponds to the quadratic variation (Da Prato and Zabczyk, 2014) of $\int_0^t \sigma(x, s) dB_s$, and it is a continuous random process of finite variation (hence time differentiable). In the particular case of a nonrandom $\sigma$, $a$ reduces then to the local variance of the noise according to the Itô isometry (Da Prato and Zabczyk, 2014). Physically, the symmetric non-negative tensor $a$ represents the (possibly random) friction coefficients of the unresolved fluid motions and the eigenvalues $\lambda_n$ have the unit of m$^2$/s.

## 2.2 Stochastic transport

The evolution law of a random tracer $\Theta$ with extensive property (e.g. temperature, salinity, buoyancy) transported by the stochastic flow, $\Theta(X_{t+\delta t}, t + \delta t) = \Theta(X_t, t)$ with $\delta t$ an infinitely small time variation, is derived by Mémin (2014); Bauer et al. (2020a) using the generalized Itô formula (Kunita, 1997). It is described by the following stochastic partial differential equation (SPDE):

$$\mathbb{D}_t \Theta = d_t \Theta + (v^* dt + \sigma dB_t) \cdot \nabla \Theta - \frac{1}{2} \nabla \cdot (a \nabla \Theta) dt = 0, \quad (2.6)$$

in which $\mathbb{D}_t$ is introduced as a stochastic transport operator and $d_t \Theta(x) := \Theta(x, t + \delta t) - \Theta(x, t)$ stands for the (forward) time-increment of $\Theta$ at a fixed point $x \in \mathcal{D}$.

This SPDE encompasses physically meaningful terms. For instance, the third term describes the tracer's advection by the unresolved flow component. As shown in Resseguier et al. (2017a); Bauer et al. (2020a), the resulting multiplicative noise $\sigma dB_t \cdot \nabla \Theta$ continuously backscatters random energy to the system through the quadratic variation $1/2 (\nabla \Theta) \cdot (a \nabla \Theta)$ of the random tracer. The last term in (2.6) depicts tracer's diffusion due to the mixing of the unresolved scales. In particular, under specific noise definitions (Mémin, 2014), the resulting diffusion can be connected to the functional eddy viscosity as introduced in many large-scale circulation models (Smagorinsky, 1963; Redi, 1982).

As an additional feature of interest, there exists an *effective* advection velocity $v^*$ in (2.6) which is defined as

$$v^* = v - \frac{1}{2} \nabla \cdot a + \sigma^*(\nabla \cdot \sigma). \quad (2.7)$$

This statistical eddy-induced velocity drift captures the action of inhomogeneity of the random field on the transported tracer and the possible divergence of the unresolved flow component. Bauer et al. (2020a) shows that the *turbophoresis* term $1/2 \nabla \cdot a$ can be interpreted as a generalization of the Stokes drift, which occurs, for example, in the Langmuir circulation (Craik and Leibovich, 1976; Leibovich, 1980).

Many useful properties of the stochastic transport operator $\mathbb{D}_t$ have been explored by Resseguier et al. (2017a, 2020); Li (2021). In particular, if a random tracer is transported by the incompressible
stochastic flow under suitable boundary conditions, then the path-wise p-th moment ($p \ge 1$) of the tracer is materially and integrally invariant, namely

$$ \mathbb{D}_t \left( \frac{1}{p} \Theta^p \right) = 0, \quad \mathrm{d}_t \int_{\mathcal{D}} \frac{1}{p} \Theta^p \mathrm{d}\mathbf{x} = 0. \qquad (2.8) $$

## 2.3 Stochastic QG model

The derivation of the stochastic geophysical models under the LU framework follows almost exactly the same path as the deterministic derivation (Vallis, 2017). In particular, a stochastic Euler-Boussinesq model has been derived by Resseguier et al. (2017b); Bauer et al. (2020a) for large-scale atmospheric and oceanic circulations, and reads:

$$ \underset{\text{Momentum equations}}{\mathbb{D}_t \mathbf{v} + f \mathbf{k} \times (\mathbf{u} \, \mathrm{d}t + \boldsymbol{\sigma}_h \, \mathrm{d}\mathbf{B}_t)} = (b \, \mathrm{d}t + \mathrm{d}b_t^\sigma) \mathbf{k} - \nabla (p \, \mathrm{d}t + \mathrm{d}p_t^\sigma), \qquad (2.9a) $$

$$ \underset{\text{Continuity equations}}{\nabla \cdot \mathbf{v}^* = 0, \quad \nabla \cdot \boldsymbol{\sigma} \, \mathrm{d}\mathbf{B}_t = 0}, \qquad (2.9b) $$

$$ \underset{\text{Thermodynamic equation}}{\mathbb{D}_t b + N^2(w^* \, \mathrm{d}t + \sigma_z \, \mathrm{d}B_t) = \frac{1}{2} \nabla \cdot (\mathbf{a}_z \, N^2) \, \mathrm{d}t}, \qquad (2.9c) $$

where $\mathbf{v} = (\mathbf{u}, w)^T$, $\boldsymbol{\sigma} \, \mathrm{d}\mathbf{B}_t = (\boldsymbol{\sigma}_h \, \mathrm{d}\mathbf{B}_t, \sigma_z \, \mathrm{d}B_t)^T$, $\mathbf{k} = (0, 0, 1)^T$, $f = f_0 + \beta y$ is the Coriolis parameter using beta-plane approximation, $p$ and $\mathrm{d}p_t^\sigma/\mathrm{d}t$ (in a distribution sense) are the time-smooth component and the fast oscillating noise of the pressure fluctuations rescaled by the background density $\rho_0$, $N^2(z) = -g \partial_z \bar{\rho}(z)/\rho_0$ is the Brunt-Väisälä frequency with $g$ the gravity value and $\bar{\rho}$ the stationary density in equilibrium that only depend on height, $b(\mathbf{x}, t) = -g \rho'(\mathbf{x}, t)/\rho_0$ is the resolved buoyancy associated with the density anomaly $\rho'$, $\mathrm{d}b_t^\sigma/\mathrm{d}t$ is a zero-mean noise that models the unresolved buoyancy fluctuations and $\mathbf{a}_z$ denotes the $z$-column vector of tensor $\mathbf{a}$. The term $\mathrm{d}b_t^\sigma$ can be seen as a thermodynamics noise coming from the nonlinear mapping law of state applied to the noise term of salinity and temperature transports. It is of the same nature as the noise term introduced in Brankart (2013).

From the system (2.9) a diverse set of approximated models under LU can be obtained through nondimensionalization and asymptotic approach with proper scaling. However, the noise introduces an additional degree of freedom that must be appropriately accounted for (Brecht et al., 2021; Bauer et al., 2020b; Resseguier et al., 2017b,c). The horizontal components of the quadratic variation are first scaled as $\mathbf{a}_h \sim \epsilon UL$, where $U$ and $L$ are typical velocity and length scales, and the factor $\epsilon$ is proportional to the ratio between the eddy kinetic energy (EKE) and the mean kinetic energy (MKE) and to the ratio between the small-scale correlation time and the large-scale advection time. From the definitions (2.2), (2.3) and (2.5), the horizontal components of the noise can be scaled as $\boldsymbol{\sigma}_h \, \mathrm{d}\mathbf{B}_t \sim \epsilon^{1/2} L$. The greater this scaling number $\epsilon$, the larger the quadratic variation, hence the stronger the path-wise noise. As shown in Resseguier et al. (2017b,c), using different levels of noise in the stochastic system allows us to model different physical regimes of the large-scale flow.

In this work, we only consider the “moderate” horizontal uncertainty, $\epsilon = O(1)$ or $\epsilon R_0 = O(R_0)$ with a small Rossby number $R_0 = U/(f_0 L) \ll 1$, in order to maintain the classical geostrophic balance.
Besides, to keep the assumption of flat isopycnal (small variation of stratification) in the QG theory, Resseguier et al. (2017b) proposed to scale the ratio between the vertical scale and the horizontal scale of the noise as $\sigma_z d\mathbf{B}_t/\sigma_h d\mathbf{B}_t \sim (H/L)(R_0/B_u)$ with the height scale $H \ll L$ and the Burger number $B_u = (L_d/L)^2 = \mathcal{O}(1)$, where $L_d$ stands for the scale of deformation radius. Moreover, the unresolved pressure and buoyancy are scaled by $dp_t^\sigma \sim \Phi = \epsilon^{1/2} f_0 L^2$ and $db_t^\sigma \sim \Phi/H$. All the resolved variables ($\mathbf{u}, w, p, b, N^2$) are scaled as in the classical framework (Vallis, 2017) and the variation of Coriolis parameter stays small, i.e. $\beta L < f_0$.

Expanding the resolved variables ($\mathbf{v}, p, b$) in power series of the small Rossby number and truncating the equations at zeroth and first orders, then restoring the dimensions of variables, the following QG system can be derived:

$$ \frac{\text{Momentum equations}}{\mathcal{D}_t^h \mathbf{u} + f_0 \mathbf{k} \times \mathbf{u}_a \, dt + \beta y \mathbf{k} \times (\mathbf{u} \, dt + \boldsymbol{\sigma}_h d\mathbf{B}_t) = -\nabla_h p_a \, dt,} \quad (2.10a) $$

$$ \frac{\text{Buoyancy equation}}{\mathcal{D}_t^h b + N^2 w_a \, dt = 0}, \quad (2.10b) $$

$$ \frac{\text{Continuity equation}}{\nabla_h \cdot \mathbf{u}_a + \partial_z w_a = 0}, \quad (2.10c) $$

where $\mathbf{u}$ and $p$ denote the geostrophic components of the resolved velocity and pressure, $(\mathbf{u}_a, w_a)^T$ and $p_a$ denote the ageostrophic components, and $\mathcal{D}_t^h$ stands for the horizontal component of the stochastic transport operator, which is defined as

$$ \mathcal{D}_t^h \Theta = d_t \Theta + (\mathbf{u}^* \, dt + \boldsymbol{\sigma}_h d\mathbf{B}_t) \cdot \nabla_h \Theta - \frac{1}{2} \nabla_h \cdot (a_h \nabla_h \Theta) \, dt \quad (2.10d) $$

with $\mathbf{u}^* := \mathbf{u} - 1/2 \nabla_h \cdot a_h$ and $\nabla_h = (\partial_x, \partial_y)^T$. The geostrophic component and the noise satisfy the following equilibrium:

$$ \frac{\text{Geostrophic balances and martingale pressure equilibrium}}{\mathbf{u} = \frac{1}{f_0} \nabla_h^\perp p, \quad \boldsymbol{\sigma}_h d\mathbf{B}_t = \frac{1}{f_0} \nabla_h^\perp dp_t^\sigma, \quad \nabla_h \cdot (\mathbf{u}^* - \mathbf{u}) = 0}, \quad (2.10e) $$

$$ \frac{\text{Hydrostatic balances}}{b = \partial_z p, \quad db_t^\sigma = \partial_z dp_t^\sigma}, \quad (2.10f) $$

where $\nabla_h^\perp = (-\partial_y, \partial_x)^T$. For sake of conciseness of the notations, we drop the horizontal subscript $h$ in the following. Combining the above equilibrium, the thermal wind balances, $\partial_z \mathbf{u} \perp \nabla_b$ and $\partial_z \boldsymbol{\sigma} d\mathbf{B}_t \perp \nabla_b db_t^\sigma$, can be deduced. In this work, we only focus on noise fully satisfying the martingale pressure equilibrium. It means that we model the mesoscale flow as a stochastic component given a large-scale flow. However, such random equilibrium is not necessary in a general setting and could be easily relaxed.

In order to represent the previous random system in terms of potential vorticity (PV) in the same way as in the classical framework (Vallis, 2017), we first take the (horizontal) curl of the first-order momentum equations (2.10a), then sum the resulting equation with the vertical derivative of
the buoyancy equation (2.10b) rescaled by $f_0/N^2$, together with the equilibrium (2.10e), (2.10f) and (2.10c), to obtain finally

### Stochastic evolution of PV

$$ \mathbb{D}_t q = \sum_{i=1,2} -J((u^* \, dt + \sigma dB_t)^i, u^i) + \frac{1}{2} \nabla \cdot (\partial_{x_i}^\perp a \nabla u^i) \, dt - \nabla \cdot (a \nabla f) \, dt \\ - \frac{f_0}{N^2} \partial_z (u^* \, dt + \sigma dB_t) \cdot \nabla b + \frac{f_0}{2N^2} \nabla \cdot (\partial_z a \nabla b) \, dt, \quad (2.11a) $$

### From PV to streamfunction

$$ q = \nabla^2 \psi + \partial_z \left( \frac{f_0}{N^2} b \right) + \beta y, \quad b = f_0 \partial_z \psi \quad (2.11b) $$

### Incompressible constraints

$$ \mathbf{u} = \nabla^{\perp} \psi, \quad \nabla \cdot \sigma d\mathbf{B}_t = \nabla \cdot (\mathbf{u}^* - \mathbf{u}) = 0, \quad (2.11c) $$

where $q$ denotes the PV, $\psi := p/f_0$ stands for the streamfunction of the resolved flow, $J(f,g) = \partial_x f \partial_y g - \partial_x g \partial_y f$ stands for the Jacobian operator and $\nabla^2 = \partial_{xx}^2 + \partial_{yy}^2$ is the two-dimensional (2D) Laplacian operator. Note that in the deterministic setting, the sources and sinks of PV in (2.11a) cancels due to the anti-symmetry of the Jacobian operator, $J(u^i, u^i) = 0$, and to the thermal wind balance, $\partial_z \mathbf{u} \perp \nabla b$. However, it should be noted that the presence of sub-grid terms in a large-scale deterministic Boussinesq system leads also in general to additional terms in the RHS of the PV equation. The additional terms in the RHS of (2.11a) are coming exactly from the same sources of difficulty inherent to large-scale flow dynamics representations (i.e. non commutation of curl operator with the sub-grid model, breaking of the Jacobian anti-symmetry as well as possibly non commutation of coarsening processes with derivatives).

## 2.4 Energy conservation and transfers

One important characteristic of the random system (2.10) or (2.11) is that it conserves path-wise (i.e. for each realization) the total energy of the resolved geostrophic flow (under natural boundary conditions), namely

$$ d_t \int_D \frac{1}{2} (|\mathbf{u}|^2 + \frac{b^2}{N^2}) d\mathbf{x} = \left( \int_D w_a b \, d\mathbf{x} - \int_D w_a b \, d\mathbf{x} \right) dt = 0. \quad (2.12) $$

The demonstration can be found in 5. We remark that this conservation property is consistent with the energy conservation of the stochastic barotropic QG system (Bauer et al., 2020a) and of the stochastic rotating shallow-water system (Brecht et al., 2021).

The energy conversions between the ensemble-mean and ensemble-eddy components are in addition explored in 5. Based on the decompositions $\mathbf{u} = \mathbb{E}[\mathbf{u}] + (\mathbf{u} - \mathbb{E}[\mathbf{u}]) := \bar{\mathbf{u}} + \mathbf{u}'$ and $b = \bar{b} + b'$, the mean kinetic energy (MKE), the eddy kinetic energy (EKE), the mean potential energy (MPE) and the eddy potential energy (EPE) are respectively defined by $1/2\langle|\bar{\mathbf{u}}|^2\rangle$, $1/2\langle|\mathbf{u}'|^2\rangle$, $1/2\langle(\bar{b}/N)^2\rangle$ and $1/2\langle(b'/N)^2\rangle$, where $\langle\cdot\rangle$ denotes the integration over domain for simplicity. The energy conversions of the proposed random model (2.10) are summarized in the following diagram.
<figure><img src="image_2.png" /></figure>

Here, $[A, B]_F = \text{Tr}(A^T B)$ stands for the Frobenius inner product of the tensors $A$ and $B$. The ensemble-mean of the path-wise energy conservation (2.12) can be recovered by summing all the conversion terms.

Compared to the energy budget of the classical deterministic system, the random system provides additional conversion terms between the mean and eddy components. In the general setting with $\sigma$ random, the classical geostrophic Reynolds stress (Bachman, 2019) is modified to $(u^*)'(u')^\dagger$ due to the possibly inhomogeneous nature of the noise. This term includes an additional term related to the Ito-Stokes drift that can be positive or negative depending on its alignment with the large-scale velocity component. The fact that $\bar{a} \, dt = \text{Var}(\sigma dB_t)$ is symmetric and non-negative definite, ensures that the integral $\langle [\nabla \bar{u}, \bar{a} \nabla \bar{u}]_F \rangle$ (resp. $\langle [\nabla \bar{b}, \bar{a} \nabla \bar{b}]_F \rangle$) is always positive, hence provides a positive transfer from MKE to EKE (resp. from MPE to EPE). When the correlation tensor $\sigma$ is nonrandom, only this term remains in addition to the term related to the classical geostrophic Reynolds stress. As the stochastic system (2.10) does not contain vertical noise component, it does not explicitly modify the conversion mechanisms between KE and PE as observed in the classical system. Nevertheless, the effect of the stochastic terms on this energy transfers is implicitly accounted via the time-integration of the resolved variables $b$ and $w_a$.

## 2.5 Girsanov transformation

The previous formulations (2.10) and (2.11) consist of only a zero-mean and temporally uncorrelated noise. However, this might not be enough in practice and including a mean or time-correlated component of the unresolved velocity field could be of crucial importance to obtain a relevant model. For instance, the eddy parameterization proposed by Zanna et al. (2017) is decomposed into a deterministic term based on the negative Laplacian of the PV transport and a zero-mean stochastic term. For the double-gyre circulation configuration, the considered deterministic parameterization allows to reproduce the eastwards jet for the coarse-resolution model, while the additional stochastic terms enhance the gyres circulation and improves the flow variability. Similarly, the random-forcing model proposed by Berloff (2005) consists in a space-time correlated stochastic process to enhance the jet extension.

The slow components of the sub-grid scales can be provided by adequate high-pass filtering of high-resolution data on the coarse grid. We aim in this work at investigating the incorporation of such
slow components within the LU framework. However, the derivation of LU models (Mémin, 2014; Resseguier et al., 2017a; Bauer et al., 2020a) relies on the martingale properties of the centered noise and we need hence to properly handle non centred Brownian terms. The Girsanov transformation (Da Prato and Zabczyk, 2014) provides a theoretical tool that fully warrants such a superposition: by a change of the probability measure, the composed noise can be centered with respect to a new probability measure while the additional drift term appears, which pulls back time-correlated sub-grid-scale components into the dynamical system. The associated mathematical description is given as follows.

Let $\Gamma_t$ be an H-valued $\mathcal{F}_t$-predictable process satisfying $\mathbb{E}_{\mathbb{P}}[\exp(\frac{1}{2} \int_0^T \|\Gamma_t\|_H^2 dt)] < +\infty$, then the process $\{\tilde{B}_t := B_t + \int_0^t \Gamma_s ds\}_{0\le t\le T}$ is an H-valued cylindrical Wiener process on $(\Omega, \mathcal{F}, \{\mathcal{F}_t\}_{0\le t\le T}, \mathbb{P})$ by a change of the probability measure:

$$ \frac{d\tilde{\mathbb{P}}}{d\mathbb{P}} = \exp\left(-\int_0^T \langle \Gamma_t, dB_t \rangle_H - \frac{1}{2} \int_0^T \|\Gamma_t\|_H^2 dt\right), \qquad (2.13) $$

where $\langle f,g \rangle_H = \int_{\mathcal{D}} (f^T g)(x) dx$ denotes the inner product on $H = (L^2(\mathcal{D}))^d$. Under this new probability measure $\tilde{\mathbb{P}}$, the previous system (2.11) must be adapted to the modified noise $\sigma d\tilde{B}_t$ with $\mathbb{E}_{\tilde{\mathbb{P}}}[\sigma d\tilde{B}_t] = 0$ by a correction of the effective drift

$$ \tilde{u}^* := u^* - \sigma \Gamma_t. \qquad (2.14) $$

Hereafter, the resulting stochastic system will be adopted with $\sigma\Gamma_t$ referred to as the Girsanov drift.

# 3 Numerical models

In this section, we first describe the discretization of the continuous stochastic QG system (2.11) (under a simplification) that will be used later for the numerical simulations. Then, we present numerical methods to estimate the spatial correlation functions of the unresolved flow component based on the spectral decomposition (2.3).

## 3.1 Discrete QG model

We consider here a vertically discretized QG model. This formulation is quite common in geophysical fluid dynamics and the derivation follows the standard methods (Vallis, 2017; McWilliams, 2006). As illustrated in Figure 1, such model consists of n ocean layers with constant thickness $H_k$ and density $\rho_k$ in each layer k. The prognostic variables $(q_k, \psi_k, u_k, \sigma_k, a_k)$ are assumed to be layer-averaged quantities.

We first apply finite differences for the buoyancy frequency $N^2$ and anomaly b on the interface between layers k and k + 1 (for k = 1, ..., n - 1):

$$ b_k = \frac{f_0(\psi_k - \psi_{k+1})}{(H_k + H_{k+1})/2}, \quad N_k^2 = \frac{g'_k}{(H_k + H_{k+1})/2}, \qquad (3.1) $$

where $g'_k := -g(\rho_k - \rho_{k+1})/\rho_0$ is introduced as the reduced gravity across the interface k. Applying subsequently finite differences for the stratification term $\partial_z(f_0 b/N^2)$ with the zero boundary conditions
at ocean surface and bottom, $b_0/N_0^2 = b_n/N_n^2 = 0$, the Helmholtz equation (2.11b) can be vertically discretized by

$$ \nabla^2 \psi_k + \frac{f_0^2}{H_k} \left( \frac{\psi_{k-1} - \psi_k}{g'_{k-1}} - \frac{\psi_k - \psi_{k+1}}{g'_k} \right) = q_k - \beta y. \qquad (3.2a) $$

In this work, we adopt the dynamical core of the QG coupled model (Q-GCM) developed by Hogg et al. (2003), which integrates numerically the state variables $q$ and $\psi$. In order to reduce the inaccuracy of numerical discretization, we do not include the PV sources and sinks (RHS terms of (2.11a)) in this work.

These terms are naturally considered when working in $u-b$ formulations such as in the primitive equations dynamical core (Bachman, 2019) or for multi-layers shallow water systems (Adcroft et al., 2019). In our case, for the type of noises considered (with moderated amplitude so as to respect the QG scaling performed in this work) these additional noise terms are small but leads to additional numerical difficulties. For stronger noise, with *de facto* other scalings with potential modification of the hydrostatic balances these additional noise terms have to be considered in the PV equation.

Taking account the numerical dissipation and the vertical entrainment $e_k$ (across interfaces), the simplified stochastic evolution of PV reads

$$ \tilde{\mathbb{D}}_t^k q_k = \left( (A_2 \nabla^2 - A_4 \nabla^4)(\nabla^2 \psi) + \frac{f_0}{H_k} (e_{k-1} - e_k) \right) dt, \qquad (3.2b) $$

where $\tilde{\mathbb{D}}_t^k$ denotes the $k$-th layer stochastic transport operator under the probability measure $\tilde{\mathbb{P}}$ (hence includes the Girsanov drift described in Section 2.5, $\nabla^4 = \nabla^2(\nabla^2)$ is the bi-Laplacian operator associated with the dissipation coefficient $A_4$.

As for the spatial discretization, we propose to use a conservative flux form $\nabla \cdot (uq)$ that corresponds exactly to the 9-points Arakawa Jacobian scheme (Arakawa and Lamb, 1981) of $J(\psi, q)$, in which the zonal and meridional components of the advection flux are defined as $2/3 (\overline{u^x} \overline{q^x}) + 1/3 (\overline{u^x} \overline{q^y})$ and $2/3 (\overline{v^y} \overline{q^y}) + 1/3 (\overline{v^x} \overline{q^d})$. Here, $\overline{\cdot}^x, \overline{\cdot}^y$ and $\overline{\cdot}^d$ denote the central averaging operations in the zonal, meridional and diagonal directions respectively. We then apply it for the advection of PV by all the additional drift (such as the noise). On the other hand, we employ simply a central winding scheme to discretize the inhomogeneous diffusion term $\nabla \cdot (a\nabla q)$. We adopt a stochastic Leapfrog scheme (Ewald and Témam, 2005) for the time-stepping of (3.2b). The inversion of the modified Helmholtz equation (2.11b) is carried out with a discrete Fourier transform method (Hogg et al., 2003).

The vertical entrainment $e$ in the classical QG models (Berloff et al., 2007; Marshall et al., 2012; Grooms et al., 2015; Zanna et al., 2017) accounts only for the Ekman pumping from the upper and bottom boundary layers, $e_0 = \nabla^\perp \cdot \tau/f_0$ and $e_n = 1/2 \delta_{ek} \nabla^2 \psi_n$, where $\tau$ denotes a given surface wind stress and $\delta_{ek}$ stands for a prescribed thickness of the bottom Ekman layer. In the Q-GCM formulation (Hogg et al., 2004), an additional entrainment $e_1$ is considered due to the Ekman pumping of mixed layer temperature. This will be further detailed in Section 4.3.

## 3.2 Parameterizations of unresolved flow

This section describes the numerical methods of the parameterization for the unresolved flow structure. We first review the empirical orthogonal functions (EOF) method which relies on an assumption of
<figure><img src="image_3.png" /><figcaption>Figure 1: Illustration for the vertical discretization of a three-layered QG system.</figcaption></figure>

quasi-stationarity of the noise covariance. To relax this stationarity constraint, we then propose a projection method enabling to update on-line the noise spatial correlations.

Let {$u_{HR}(x, t_i)}_{i=1,...,N}$ be the finite set of velocity snapshots provided by a high-resolution (HR) simulation. We first build the spatial local fluctuations $u_f(x, t_i)$ of each snapshot at the coarse-grid points. In particular, for the QG system, one can first perform a high-pass filtering with a 2D Gaussian convolution kernel $G$ on each HR streamfunction $\psi_{HR}$, to obtain the streamfunction fluctuations, $\psi_f(x, t_i) = ((I - G) \star \psi_{HR})(x, t_i)$ (only for the coarse-grid points $x$). Then, the geostrophic velocity fluctuations can be derived by $u_f = \nabla_{\perp} \psi_f$.

We next center the data set by $u'_f = u_f - \overline{u_f}^t$ (with $\overline{\cdot}^t$ the temporal mean) and perform the EOF procedure (Sirovich, 1987) to get a set of orthogonal temporal modes {$\alpha_n$}$_{n=1,...,N}$ and orthonormal spatial modes {$\phi_n$}$_{n=1,...,N}$ satisfying

$$u'_{f}(\mathbf{x}, t_i) = \sum_{n=1}^{N} \alpha_n(t_i) \phi_n(\mathbf{x}), \qquad (3.3)$$

where $\langle \phi_m, \phi_n \rangle_H = \delta_{m,n}$, $\overline{\alpha_m \alpha_n}^t = \lambda_n \delta_{m,n}$, $H := L^2(D)$ and $D = A \times [-H, 0]$ denotes the 3D domain with the horizontal area $A$ and the total depth $H$.

Truncating subsequently the modes (with $M \ll N$) and rescaling them by a small-scale decorrelation time $\tau$, the $k$-th layer noise and its variance are built by

$$\sigma_k(\mathbf{x}) \mathrm{d}\tilde{B}_t = \tau^{1/2} \sum_{n=1}^{M} \lambda_n^{1/2} \phi_{n,k}(\mathbf{x}) \mathrm{d}\beta_n(t), \qquad (3.4a)$$

$$a_k(\mathbf{x}) = \tau \sum_{n=1}^{M} \lambda_n \phi_{n,k}(\mathbf{x}) \phi_{n,k}^{T}(\mathbf{x}). \qquad (3.4b)$$

Note that this time scale $\tau$ is used to match the fact that the noise in (2.1) admits unit of meter. In practice, we often consider the coarse-grid simulation time-step $\Delta t_{LR}$.
In this work, we simply propose to represent the Girsanov drift by the temporal mean of the unresolved eddies, that is

$$ \sigma_k(\mathbf{x})\Gamma_t = \sum_{n=1}^{N} \langle \overline{u}_f^t, \phi_n \rangle_{\mathcal{H}} \phi_{n,k}(\mathbf{x}). \qquad (3.4c) $$

We remark that in future works this stationary drift term could be generalized to slowly varying fluctuation components or to drift correction with respect to observation (Dufée et al., 2022).

The previous EOF procedure is an efficient off-line learning method, yet it relies on a strong stationary assumption, and hence leads to a sequence of random flow fields with no connection with the resolved dynamics. In the following, we propose an novel approach that project the EOF-based noise on a resolved state variable. In order to enforce the noise to act only on the resolved momentum (2.10a) and provide efficient backscattering of KE, we consider constraining the noise along the iso-surfaces of the stratification, $\tilde{\sigma}d\tilde{B}_t \perp \nabla\partial_z b$ with $\tilde{b} := f_0b/N^2$. This projection procedure reads

$$ \tilde{\sigma}_k(\mathbf{x}, t) d\tilde{B}_t = A^{1/2} \mathbf{P}_k(\mathbf{x}, t) \sigma_k(\mathbf{x}) d\tilde{B}_t, \quad \mathbf{P}_k := \mathbf{I}_2 - \frac{\nabla(\partial_z \tilde{b})_k (\nabla(\partial_z \tilde{b})_k)^T}{|\nabla(\partial_z \tilde{b})_k|^2}, \qquad (3.5) $$

where $\sigma_k d\tilde{B}_t$ is generated as in (3.4a) and A is a scaling factor to ensure that the noise amplitude is projection invariant. From the definition of the noise quadratic variation (2.5), tensor a defined through (3.4b) is updated as $\tilde{a}_k = A \mathbf{P}_k a_k \mathbf{P}_k^T$. In the next section, we show that the resulting non-stationary noise leads to significant improvement of variability for coarse-resolution simulations.

# 4 Numerical results

In this section, we discuss and compare numerical results of several models run within different configurations. The objective consists in improving the variability of low-order models at low Reynolds numbers flow, and in the particular case of unresolved baroclinic instabilities.

## 4.1 Model configurations

We consider here a finite box ocean at mid-latitude driven by an idealized (stationary and symmetric) wind stress $\tau = (-\tau_0 \cos(2\pi y)/L_y, 0)^T$. A mixed horizontal boundary condition is used for the k-th layer streamfunction: $\psi_k|_{\partial A} = f_k(t)$ and $\partial_n^2 \psi_k|_{\partial A} = -(\alpha_{bc}/\Delta x)\partial_n \psi_k|_{\partial A}$ (same for the 4-th order derivative). Here, $f_k$ is a time-dependent function constrained by mass conservation (McWilliams, 1977), $\Delta x$ stands for the horizontal resolution and $\alpha_{bc}$ is a nondimensional coefficient associated to the slip boundary conditions (Haidvogel et al., 1992).

A deterministic eddy-resolving ($\Delta x = 5$ km) model is first simulated and referred to as the references (REF). It is then compared to several coarse-resolution ($\Delta x = 40$ km, 80 km, 120 km) models: the benchmark deterministic model (DET), two stochastic models with a stationary EOF-based noise (STO-EOF) and a non-stationary projection noise (STO-EOF-P). The eddy-resolving model starts from a quiescent initial condition, whereas spin-up conditions downsampled from REF (after 90-years integration) are adopted for all the coarse-resolution models. The common parameters for all the simulations are listed in Table 1, whereas resolution dependent parameters are presented separately in Table 2.
<table><caption>Table 1: Common parameters for all the models.</caption><thead><tr><th>Parameters</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td>X × Y</td><td>(3840 × 4800) km</td><td>Domain size</td></tr><tr><td>H<sub>k</sub></td><td>(350, 750, 2900) m</td><td>Mean layer thickness</td></tr><tr><td>ρ</td><td>1000 kg m<sup>3</sup></td><td>Density</td></tr><tr><td>g'<sub>k</sub></td><td>(0.025, 0.0125) m s<sup>-2</sup></td><td>Reduced gravity</td></tr><tr><td>δ<sub>ek</sub></td><td>2 m</td><td>Bottom Ekman layer thickness</td></tr><tr><td>τ<sub>0</sub></td><td>2 × 10<sup>-5</sup> m<sup>2</sup> s<sup>-2</sup></td><td>Wind stress magnitude</td></tr><tr><td>α<sub>bc</sub></td><td>0.2</td><td>Mixed boundary condition coefficient</td></tr><tr><td>f<sub>0</sub></td><td>9.375 × 10<sup>-5</sup> s<sup>-1</sup></td><td>Mean Coriolis parameter</td></tr><tr><td>β</td><td>1.754 × 10<sup>-11</sup> (m s)<sup>-1</sup></td><td>Coriolis parameter gradient</td></tr><tr><td>L<sub>d</sub></td><td>(39, 22) km</td><td>Baroclinic Rossby radii</td></tr></tbody></table>

<table><caption>Table 2: Values of grid varying parameters.</caption><thead><tr><th>Δx (km)</th><th>Δt (s)</th><th>A<sub>4</sub> (m<sup>4</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>5</td><td>600</td><td>2 × 10<sup>9</sup></td></tr><tr><td>40</td><td>1200</td><td>5 × 10<sup>11</sup></td></tr><tr><td>80</td><td>1440</td><td>5 × 10<sup>12</sup></td></tr><tr><td>120</td><td>1800</td><td>1 × 10<sup>13</sup></td></tr></tbody></table>

Snapshots of the relative vorticity on the two top layers are shown in Figure 2. The dynamics of REF model is mainly characterized by the energetic eastward jet with adjacent recirculations and by the fast westward Rossby waves. The former results from the most active mesoscale eddies effect through baroclinic instability. However, this instability cannot be resolved once the horizontal resolution is similar or greater than the baroclinic deformation radius maximum (39 km here). For instance, without any eddy parameterization, the DET (40 km) simulation only generates a smooth field due to an over-dissipation effect.

To calibrate the structure of the unresolved flow (noise and Girsanov drift), the REF data are sampled at each day during the first ten years after the spin-up, this is followed by the procedure presented in Section 3.2 with a fixed 200 km wide Gaussian filter. Considering the computational cost, we employ only 50 EOFs to build the noises (3.4a) for the stochastic simulations at different resolutions. However, the eigenvalues of these truncated EOFs are amplified in order to capture 95 % energy of the full set of modes. An example of the first EOF ($\sqrt{\lambda_1}\phi_1$) and the noise variance ($a$ defined in (3.4b)) is illustrated by Figure 3. The noise structure is mainly characterized by the most energetic fluctuations in the jet region.

Including the additional advection of PV by the calibrated unresolved flow component, both STO-EOF and STO-EOF-P models are able to reproduce the eastward jet on the coarse mesh, and the latter seems to enhance fluctuations everywhere. For instance, such result is illustrated at a resolution of 40 km by Figure 2. The ability of these coarse models to reproduce the statistical properties of the REF model will be diagnosed and analyzed more precisely in the next section.
<figure><img src="image_4.png" /><figcaption>Figure 2: Snapshots of relative vorticity (divided by $f_0$) for the two upper layers (by rows) provided by different simulations (by columns) after 60-years integration.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Illustration for noise's first EOF (zonal and meridional components in first two sub-figures) and for noise's variance (diagonal components in last two sub-figures) at resolution 40 km.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Comparison of time-mean contour of barotropic (top row) and first baroclinic (bottom row) streamfunctions for different models (by columns). Note that 1 Sv = $10^6 \text{ m}^3\text{/s}$.</figcaption></figure>

## 4.2 Statistical diagnostics

To investigate the prediction of the long-term statistics, we sample the data (streamfunction) of each run at 15-day intervals over the 120 years (after the spin-up) and compute the statistics over the last 100 years.

We first compare the temporal mean of the modal streamfunctions for both coarse models to that subsampled (in space) of the eddy-resolving model. From Figure 4, we observe that both barotropic and baroclinic modal mean of the *REF* characterize the eastwards jet, and the two stochastic models (at 40 km) enable to reproduce qualitatively the local structures of both vertical modes predicted by the *REF*. Conversely, the *DET* model can only capture the symmetric double-gyre structure.

To quantify the models’ statistics at different resolutions, we then employ the spatial root-mean-square error (RMSE) of the temporal mean and the standard deviation (std) between each coarse model and the subsampled *REF* one. To assess in a single measure the mean and the variance reconstruction, we adopt the Gaussian relative entropy (GRE) proposed by Grooms et al. (2015) that
<figure><img src="image_7.png" /><figcaption>Figure 5: Comparison of variability measures (by sub-figures) for different coarse models (by colors and groups in each sub-figure). Note that the y-axis of the last sub-figure is in log-scale due to poor variability of the DET models.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Temporal mean of vertically integrated KE spectra for different models (by colors) at different resolutions (by sub-figures).</figcaption></figure>

is defined as

$$ \text{GRE} = \frac{1}{|\mathcal{D}|} \int_{\mathcal{D}} \frac{1}{2} \left( \frac{(\overline{\psi^M}^t - \overline{\psi^R}^t)^2}{\sigma^2[\psi^M]} + \frac{\sigma^2[\psi^R]}{\sigma^2[\psi^M]} - 1 - \log \left( \frac{\sigma^2[\psi^R]}{\sigma^2[\psi^M]} \right) \right) dx, \quad (4.1) $$

where $\sigma$ denotes the std, $\psi^M$ (resp. $\psi^R$) stands for the streamfunction of the coarse model (resp. for that subsampled from the REF).

These three global criterion are computed for all the coarse models at different resolutions and the results are summarized in Figure 5. For all the resolutions tested the two random models have lower errors than the DET model both in terms of mean and std. In addition, the STO-EOF-P models enable to reduce the RMSE of the std and the GRE, which means an improvement on the prediction of variability.

We investigate subsequently the models’ variability based on the energy analysis. The time-averaged KE spectra provided by different models are compared in Figure 6. Unsurprisingly, compared to the REF, the DET coarse models produce a severe lack of resolved KE (Arbic et al., 2013; Kjellsson and Zanna, 2017) due to the excessive dissipation without any eddy parameterization. Both stochastic models bring more energy over all the wavenumbers, particularly at large scales. Moreover, when the resolution increases (from 120 km to 40 km), their spectra slopes in the inertial range become more and more closer to that of the REF model.

We compare finally the energy decomposition as described in Section 2 for all the models. However, we focus only on a single realization here and the eddy-mean decomposition of $u$ and $b$ is simply build
<figure><img src="image_9.png" /><figcaption>Figure 7: Comparison of the global MKE, EKE, MPE and EPE (by sub-figures) for different coarse models (by colors and groups in each sub-figure). The corresponding values of the REF model: MKE = 0.23 J/m³, EKE = 0.62 J/m³, MPE = 4.01 J/m³ and EPE = 0.34 J/m³.</figcaption></figure>

from a global time-average. From Figure 7, we observe that both random models have higher MKE and MPE than the DET model. In particular, at resolution 40 km, they are almost at the same order as that of the REF model. Compared to the STO-EOF models, STO-EOF-P models produce globally higher EKE and EPE at all the different resolutions. However, they remain much lower than that of the REF ones.

Furthermore, the horizontal density of the EKE and EPE provided by different models are compared in Figure 8. It shows that the STO-EOF-P model improves locally the eddy energy almost everywhere except in the jet region. A more precise parameterization of the unresolved flow is required to explore in future works to improve locally the energy transfers.

In summary, the diagnostic results of this section demonstrate that including the unresolved flow structure inferred from the data into the proposed random models enable to reproduce on the coarse mesh the mean flow of the eddy-resolving simulation. In addition, as performed here through a projection onto the iso-surfaces of vertical stratification, providing an adequate dynamics for the unresolved flow structure is of crucial importance to increase the variability of the coarse models.

## 4.3 Including mixed layer dynamics

In this section, we follow further the Q-GCM formulation (Hogg et al., 2004) to embed a mixed layer (of thickness $H_m = 100$ m) in the upper QG layer. This brings an additional diabatic forcing into the PV equation (3.2b) due to the Ekman pumping of sea surface temperature (SST). The vertical velocity across the first ocean interface is $e_1 = -1/2 (\Delta_m T / \Delta_1 T) e_0$, where $\Delta_m T = T_m - T_1$ (resp. $\Delta_1 T = T_1 - T_2$) denotes the temperature difference across the bottom boundary of the mixed layer (resp. of the upper QG layer). Here, the QG potential temperature are fixed to be (287, 282, 276) K.

The additional forcing $e_1$ evolves in time as the SST does. In the present work, we only considered the geostrophic component of the unresolved noise in the mixed layer. In this case, the stochastic evolution of SST reads

$$ \tilde{\mathbb{D}}_t^m T_m + T_m \nabla \cdot \mathbf{u}_m \, dt = \left( (K_2 \nabla^2 - K_4 \nabla^4) T_m + \frac{T_m + T_1}{2H_m} e_0 \right) \, dt, \quad (4.2) $$

where $\mathbf{u}_m = \nabla^\perp \psi_1 - \tau^\perp / (f_0 H_m)$ is the divergent mixed layer velocity, $\tilde{\mathbb{D}}_t^m$ is the mixed layer stochastic
<figure><img src="image_10.png" /></figure>

Figure 8: Comparison of the vertically-integrated EKE density (top row) and EPE density (bottom row) for different models (by columns). Note that the colorbar of coarse-resolution (40 km) models is different than that of the *REF* model (5 km).
<figure><img src="image_11.png" /><figcaption>Figure 9: Snapshots of relative vorticity (divided by $f_0$) for the two upper layers (by rows) provided by different simulations (by columns) after 60-years integration.</figcaption></figure>

transport operator under probability measure $\tilde{\mathbb{P}}$ (associated to $u_m$), $K_2 = 100 \text{ m}^2/\text{s}$ and $K_4 = A_4$ (see values in Table 2). Note that this shallow-water type equation is derived by vertical integration of a 3D heat equation. Equation (4.2) is numerically discretized by the central winding scheme (in space) and the stochastic Leapfrog scheme (in time).

As shown in Figure 9, the inclusion of the additional forcing coming from the SST evolution generates more small-scale eddies almost everywhere in REF. The improvements of the two stochastic models observed in the previous test case still hold here, namely the zonal jet is reproduced and perturbed. In addition, the variability of the gyres seems also to be enhanced in this case.

Subsequently, an analogue procedure is performed for all the coarse-resolution models in this configuration. For instance, the results of the time-mean streamfunctions at 40 km are provided in Figures 10 and a conclusion consistent with the previous test case is recovered.

As illustrated in Figure 11, the STO-EOF-P models, with a noise along the iso-surfaces of vertical stratification, involves an additional SST random forcing that brings an higher KE backscattering than the stationary noise model at both resolutions. This seems to highlight the importance of the non-stationary characteristic of the noise brought by the projection.
<figure><img src="image_12.png" /><figcaption>Figure 10: Comparison of time-mean contour of barotropic (top row) and first baroclinic (bottom row) streamfunctions for different models (by columns).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 11: Temporal mean of vertically integrated KE spectra for different models (by colors) at different resolutions (by sub-figures).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: Comparison of variability measures (by sub-figures) for different coarse models (by colors and groups in each sub-figure). The EKE of the *REF* model is 8.26 J/m³.</figcaption></figure>

As shown in Figure 12, the *STO-EOF-P* models provide the best results for all the metrics of variability at different resolutions. The local structures of the eddy energy are also better than those recovered in the other coarse models, which is demonstrated by Figure 13.

# 5 Conclusions

A stochastic parameterization of the unresolved eddy field is successfully implemented in a well-known QG model. The empirical spatial correlation of the small-scale noise has been first estimated from eddy-resolving simulation data. An additional correction drift that is fully justified in the stochastic setting considered and ensuing from a change of probability measure has been introduced. This non intuitive term seems quite important in the reproduction of the eastward jet within the wind-driven double-gyre circulation. In addition, a projection method has been proposed to constrain the noise living along the iso-surfaces of the vertical stratification. The resulting non-stationary noise model enables us to improve the intrinsic variability of the large-scale resolved flow. This improvement has been demonstrated through some statistical criterion. We have highlighted that both stationary and non-stationary noise models provide improvements at ocean climatic scale in terms of variability metrics, compared to a corresponding deterministic coarse model. At that resolution, the effects of the mesoscale eddies within the ocean basin at midlatitude (of deformation radius around 30 ~ 40 km) are better represented by the random models, even though the baroclinic instability can not be resolved.

The numerical results presented in this work encourage us to implement the proposed random model on more complex and realistic flow configurations. A following study for a fully coupled model of the ocean and atmosphere (Hogg et al., 2003) is already in progress. Including random forcing into the bottom layer of the atmosphere leads to an unsteady wind stress and to a random diabatic forcing due to Ekman pumping of the sea surface temperature (SST) in both atmospheric and oceanic mixed layers. Thus, it will be interesting to verify the response of the ocean dynamics in terms of variability and energy to these highly variable forcing, and particularly to compare the results to a stationary wind forced ocean dynamics as well as to an unforced dynamics. The stochastic primitive Boussinesq models (2.9) could then be tested. However, in that case the random pressure term in the LU model has to be a priori modeled or computed based on proper scaling of hydrostatics.

Another possible research axis focuses on new parameterization methods for the unresolved flow
<figure><img src="image_15.png" /><figcaption>Figure 13: Comparison of the vertically-integrated EKE density (top row) and EPE density (bottom row) for different models (by columns).</figcaption></figure>
component. On the one hand, techniques relying on the dynamic mode decomposition (Schmid, 2010; Gugole and Franzke, 2020), or on the spectral analysis of the Koopman operator (Mezić, 2013; Giannakis, 2019) could be explored in practice to estimate the dynamics of the noise and the Girsanov drift from high-resolution data. The objective will consist in evaluating their ability to characterize the long-terms impacts of some small-scale events.

Preliminary very encouraging results for the same stochastic QG model as the one studied in the present work can be found in Li et al. (2022). On the other hand, it is also interesting to design adequate noises that provide efficient energy conversion (Bachman, 2019; Jansen et al., 2019) based on the analysis presented in this work. The objective is to release the APE through the noise quadratic variation (2.10b) and to backscatter the unresolved KE to the resolved KE through the advection of momentum by the noise (2.10a).

## Data Availability Statement

The source codes to reproduce the results are available at https://github.com/matlong/qgcm_lu, all the coarse-resolution data can be found in https://doi.org/10.5281/zenodo.6818751, and the high-resolution data can be either reproduced by the restart files or accessed by contacting the corresponding author.

## Acknowledgments

The authors acknowledge the support of the ERC EU project 856408-STUOD. We would like to thank the technical support provided by Charles Deltel and Pranav Chandramouli.

## Appendix A. Conservation of path-wise energy

In this section, we show that the stochastic Boussinesq QG system (2.9), without considering the beta effect on the noise, preserves along time the global energy of the resolved geostrophic flow component for each realization. To this end, we first recall the Itô's integration-by-part formula Kunita (1997): let $d_t X = F dt + G dB_t$ and $d_t Y = F' dt + G' dB_t$ be two S(P)DEs driven by a Brownian motion B satisfying $E[GdB_t] = E[G'dB_t] = 0$, then

$$ 
\begin{aligned}
d_t(XY) &= Xd_tY + Yd_tX + d\langle X, Y \rangle_t \\
&= Xd_tY + Yd_tX + GG' dt,
\end{aligned} 
\quad (5.1) $$

where $\langle X, Y \rangle_t = \lim_{n \to +\infty} \sum_{i=1}^{p_n} (X_{t_i}^n - X_{t_{i-1}}^n)(Y_{t_i}^n - Y_{t_{i-1}}^n)$ stands for the joint quadratic variation between the processes X and Y with $0 = t_0^n < t_1^n < \dots < t_{p_n}^n = t$ a partition of the interval $[0, t]$, and the limit, if it exists, is defined in the sense of convergence in probability. In particular, we have $\langle B, B \rangle_t = t$.

We now expand equation (2.10a) for the resolved geostrophic velocity $u = (u, v)^T$ and equation
(2.10b) for the resolved buoyancy b,

$$d_t u = -(\boldsymbol{u}^* \, dt + \boldsymbol{\sigma} \mathrm{d}B_t) \cdot \boldsymbol{\nabla} u + \left(\frac{1}{2} \boldsymbol{\nabla} \cdot (\boldsymbol{a} \boldsymbol{\nabla} u) + f_0 v_a + \beta y v - \partial_x p_a\right) dt \quad (5.2a)$$

$$d_t v = -(\boldsymbol{u}^* \, dt + \boldsymbol{\sigma} \mathrm{d}B_t) \cdot \boldsymbol{\nabla} v + \left(\frac{1}{2} \boldsymbol{\nabla} \cdot (\boldsymbol{a} \boldsymbol{\nabla} v) - f_0 u_a - \beta y u - \partial_y p_a\right) dt \quad (5.2b)$$

$$d_t b = -(\boldsymbol{u}^* \, dt + \boldsymbol{\sigma} \mathrm{d}B_t) \cdot \boldsymbol{\nabla} b + \frac{1}{2} \boldsymbol{\nabla} \cdot (\boldsymbol{a} \boldsymbol{\nabla} b) \, dt - N^2 w_a \, dt, \quad (5.2c)$$

together with the geostrophic equilibrium (2.10e), the hydrostatic balances (2.10f) and the continuity equation (2.10c). Note that $\boldsymbol{a} \, dt$ corresponds to the quadratic variation of the noise $\boldsymbol{\sigma} \mathrm{d}B_t$, which is $\boldsymbol{a}$ defined in (2.5). Applying equation (5.1) for the density of APE, $b^2/(2N^2)$, and integrating subsequently over the ocean domain $\mathcal{D} = \mathcal{A} \times [0, -H]$, we have

$$\begin{aligned}
\int_{\mathcal{D}} d_t \left(\frac{b^2}{2N^2}\right) d\boldsymbol{x} &= \int_{\mathcal{D}} \frac{1}{N^2} \left(b d_t b + \frac{1}{2} d\langle b, b \rangle_t\right) d\boldsymbol{x} \\
&= \int_{\mathcal{D}} \frac{1}{N^2} \left( -b(\boldsymbol{u}^* \, dt + \boldsymbol{\sigma} \mathrm{d}B_t) \cdot \boldsymbol{\nabla} b + \frac{1}{2} b \boldsymbol{\nabla} \cdot (\boldsymbol{a} \boldsymbol{\nabla} b) \, dt \right. \\
&\qquad \left. - N^2 w_a b \, dt + \frac{1}{2} (\boldsymbol{\nabla} b) \cdot \boldsymbol{a} \boldsymbol{\nabla} b \, dt \right) d\boldsymbol{x} \\
&= - \int_{\mathcal{D}} \boldsymbol{\nabla} \cdot \left( (\boldsymbol{u}^* \, dt + \boldsymbol{\sigma} \mathrm{d}B_t) \frac{b^2}{2N^2} \right) d\boldsymbol{x} - \int_{\mathcal{D}} w_a b \, d\boldsymbol{x} \, dt \\
&= - \int_{\mathcal{D}} w_a b \, d\boldsymbol{x} \, dt,
\end{aligned} \quad (5.3)$$

where the third equality results from the geostrophic balances (2.10e) and the last equality comes from the divergence theorem with ideal boundary conditions, $\boldsymbol{u} \cdot \boldsymbol{n} |_{\partial \mathcal{A}} = \boldsymbol{\sigma} \mathrm{d}B_t \cdot \boldsymbol{n} |_{\partial \mathcal{A}} = 0$ (which leads to $\boldsymbol{u}^* \cdot \boldsymbol{n} |_{\partial \mathcal{A}} = 0$ from equations (2.7) and (2.5)). Similarly to the previous calculation, one can show that

$$\int_{\mathcal{D}} d_t \left(\frac{u^2}{2}\right) d\boldsymbol{x} = \int_{\mathcal{D}} u(f_0 v_a + \beta y v - \partial_x p_a) d\boldsymbol{x} \, dt, \quad (5.4a)$$

$$\int_{\mathcal{D}} d_t \left(\frac{v^2}{2}\right) d\boldsymbol{x} = - \int_{\mathcal{D}} v(f_0 u_a + \beta y u + \partial_y p_a) d\boldsymbol{x} \, dt. \quad (5.4b)$$

Summing these two equations, we deduce the evolution of KE:

$$\begin{aligned}
\int_{\mathcal{D}} d_t \left(\frac{|\boldsymbol{u}|^2}{2}\right) d\boldsymbol{x} &= - \int_{\mathcal{D}} \left( (\boldsymbol{\nabla}^{\perp} p) \cdot \boldsymbol{u}_a^{\perp} + \boldsymbol{\nabla} \cdot (\boldsymbol{u} p_a) \right) d\boldsymbol{x} \, dt \\
&= - \int_{\mathcal{D}} \left( \boldsymbol{\nabla}^{\perp} \cdot (p \boldsymbol{u}_a^{\perp}) - p \boldsymbol{\nabla} \cdot \boldsymbol{u}_a \right) d\boldsymbol{x} \, dt \\
&= - \int_{\mathcal{A}} \int_{-H}^{0} p \partial_z w_a \, dz \, d\mathcal{A} \, dt \\
&= - \int_{\mathcal{A}} \left( [p w_a]_{-H}^{0} - \int_{-H}^{0} w_a \partial_z p \, dz \right) d\mathcal{A} \, dt \\
&= \int_{\mathcal{D}} w_a b \, d\boldsymbol{x} \, dt,
\end{aligned} \quad (5.5)$$
where the first equality comes from the geostrophic balances (2.10e), the third equality results from both the Stokes theorem with ideal boundary condition of $u_a$ and the continuity equation (2.10c), the fourth equality derives from the vertical integration-by-parts formula, and the last equality is based on both the vertical boundary conditions, $w_a(z = 0) = w_a(z = -H) = 0$ and the hydrostatic balances (2.10f). Finally, by summing equations (5.5) and (5.3), one show that the stochastic system (5.2) conserves the path-wise total energy of the resolved geostrophic flow component, namely

$$d_t \int_D \frac{1}{2} (|\mathbf{u}|^2 + \frac{b^2}{N^2}) d\mathbf{x} = \int_D d_t (\frac{|\mathbf{u}|^2}{2}) d\mathbf{x} + \int_D d_t (\frac{b^2}{2N^2}) d\mathbf{x} = 0. \quad (5.6)$$

Since this conservation property is path-wise, it holds also for their ensemble-mean.

## Appendix B. Conversion of ensemble energies

In this section, we show briefly the conversions between the ensemble-variance and the energy of the ensemble-mean. To that end, we first decompose each prognostic variable $\Theta$ (could be $u, v, b$ or $a$) of (5.2) into a mean component $\bar{\Theta} := \mathbb{E}[\Theta]$ and an eddy component $\Theta' := \Theta - \mathbb{E}[\Theta]$. The other time-smooth components $F dt$ apart from the geostrophic transport operator are decomposed in the same way. For instance, $F = -N^2(\overline{w_a} + w_a')$ in the buoyancy equation. However, the noise is only an eddy component since $\mathbb{E}[\sigma d\mathbf{B}_t] = 0$. Substituting such decomposition into equations (5.2) and taking subsequently the ensemble-mean of the expanded equation, we obtain

$$\partial_t \bar{\Theta} = \bar{F} - \overline{u^*} \cdot \nabla \bar{\Theta} + \frac{1}{2} \nabla \cdot (\overline{a} \nabla \bar{\Theta}) - (\overline{u^*})' \cdot \nabla \bar{\Theta}' + \frac{1}{2} \nabla \cdot (\overline{a'} \nabla \bar{\Theta}') \quad (5.7a)$$

where the second and the third terms on the RHS describe the mean-mean interactions, whereas the last two terms represent the effect of eddy-eddy interactions on the ensemble-mean component. Under natural boundary conditions, one can easily show the time evolution of the global energy for $\bar{\Theta}$:

$$\frac{d}{dt} \int_D \frac{\bar{\Theta}^2}{2} d\mathbf{x} = \int_D \bar{F} \bar{\Theta} d\mathbf{x} + \int_D \overline{(u^*)' \Theta'} \cdot \nabla \bar{\Theta} d\mathbf{x} \\ - \frac{1}{2} \int_D (\nabla \bar{\Theta}) \cdot (\overline{a} \nabla \bar{\Theta} + \overline{a'} \nabla \bar{\Theta}') d\mathbf{x}. \quad (5.7b)$$

Applying this result for $\bar{u}$ and $\bar{v}$ (resp. for $\bar{b}$) together with (5.5) (resp. with (5.3)), one can deduce the evolution of the MKE (resp. MPE), and the conversion terms are summarized in the diagram in Section 2.4.

Subtracting next (5.7a) from (5.2), we deduce the evolution of $\Theta'$:

$$\begin{aligned} d_t \Theta' = & - (\overline{u^*} \cdot \nabla \Theta' - \frac{1}{2} \nabla \cdot (\overline{a} \nabla \Theta')) dt \\ & - ((u^*)' dt + \sigma d\mathbf{B}_t) \cdot \nabla \bar{\Theta} - \frac{1}{2} \nabla \cdot (\overline{a'} \nabla \bar{\Theta}) dt + F' dt \\ & - ((u^*)' dt + \sigma d\mathbf{B}_t) \cdot \nabla \Theta' - \frac{1}{2} \nabla \cdot (\overline{a'} \nabla \Theta') dt \\ & + (\overline{(u^*)'} \cdot \nabla \Theta' - \frac{1}{2} \nabla \cdot (\overline{a'} \nabla \Theta')) dt, \end{aligned} \quad (5.8a)$$
where the first two lines describe the eddy-mean interactions, whereas the last two lines represent the effect of eddy-eddy interactions on the ensemble-eddy component. To evaluate the ensemble-mean of the energy for $\Theta'$, the Itô's integration-by-part formula (5.1) is required again. In this case, one can first show that the mean of the joint quadratic variation process for $\Theta'$ reads

$$ \overline{d\langle\Theta', \Theta'\rangle_t} = (\overline{\nabla\Theta'}) \cdot (\overline{a\nabla\Theta'}) + (\overline{\nabla\Theta'}) \cdot (\overline{a'\nabla\Theta'}) + (\overline{\nabla\Theta'}) \cdot (\overline{a'\nabla\Theta'}) \\ + (\overline{\nabla\Theta'}) \cdot (\overline{a'\nabla\Theta'}) + (\overline{\nabla\Theta'}) \cdot (\overline{a\nabla\Theta}). \qquad (5.8b) $$

Subsequently, we can deduce the evolution of the global variance for $\Theta$, namely

$$ \frac{d}{dt} \int_D \frac{(\Theta')^2}{2} dx = \int_D \overline{F'\Theta'} dx - \int_D \overline{(u^*)'\Theta'} \cdot \nabla \overline{\Theta} dx \\ + \frac{1}{2} \int_D \nabla(\overline{\Theta}) \cdot (\overline{a\nabla\Theta} + \overline{a'\nabla\Theta'}) dx. \qquad (5.8c) $$

Applying similarly this result for $\bar{u}$ and $\bar{v}$ (resp. for $\bar{b}$), we can deduce the evolution of the EKE (resp. EPE), and the results are shown in the diagram in Section 2.4.

## References

Adcroft, A., Anderson, W., Balaji, V., Blanton, C., Bushuk, M., Dufour, C.O., Dunne, J.P., Griffies, S.M., Hallberg, R., Harrison, M.J., Held, I.M., Jansen, M.F., John, J.G., Krasting, J.P., Langenhorst, A.R., Legg, S., Liang, Z., McHugh, C., Radhakrishnan, A., Reichl, B.G., Rosati, T., Samuels, B.L., Shao, A., Stouffer, R., Winton, M., Wittenberg, A.T., Xiang, B., Zadeh, N., Zhang, R., 2019. The GFDL global ocean and sea ice model OM4.0: Model description and simulation features. Journal of Advances in Modeling Earth Systems 11, 3167–3211.

Arakawa, A., Lamb, V.R., 1981. A potential entrophy and energy conserving scheme for the shallow water equations. Monthly Weather Review 109, 18–36.

Arbic, B.K., Polzin, K.L., Scott, R.B., Richman, J.G., Shriver, J.F., 2013. On eddy viscosity, energy cascades, and the horizontal resolution of gridded stallite altimeter products. Journal of Physical Oceanography 43, 283–300.

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

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

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

Bauer, W., Chandramouli, P., Li, L., Mémin, E., 2020b. Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models. Ocean Modelling 151, 101646.
Berloff, P., 2005. Random-forcing model of the mesoscale oceanic eddies. Journal of Fluid Mechanics 529, 71–95.

Berloff, P., Hogg, A.M., Dewar, W.K., 2007. The turbulent ocsillator: A mechanism of low-frequency variability of the wind-driven ocean gyres. Journal of Physical Oceanography 37, 2363–2386.

Brankart, J., 2013. Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modeling. Ocean Modelling 66, 64–76.

Brecht, R., Li, L., Bauer, W., Mémin, E., 2021. Rotating shallow water flow under location uncertainty with a structure-preserving discretization. Journal of Advances in Modeling Earth Systems 13.

Cotter, C., Crisan, D., Holm, D.D., Pan, W., Shevchenko, I., 2019. Numerically modeling stochastic Lie transport in fluid dynamics. Multiscale Modeling & Simulation 17, 192–232.

Cotter, C., Crisan, D., Holm, D.D., Pan, W., Shevchenko, I., 2020. Modelling uncertainty using stochastic transport noise in a 2-layer quasi-geostrophic model. Foundations of Data Science 2, 173.

Craik, A., Leibovich, S., 1976. A rational model for Langmuir circulations. Journal of Fluid Mechanics 73, 401–426.

Da Prato, G., Zabczyk, J., 2014. Stochastic equations in infinite dimensions. Encyclopedia of Mathematics and its Applications. 2 ed., Cambridge University Press.

Dufée, B., Mémin, E., Crisan, D., 2022. Stochastic parametrization: An alternative to inflation in ensemble Kalman filters. Quarterly Journal of the Royal Meteorological Society 148, 1075–1091.

Ewald, B., Témam, R., 2005. Numerical analysis of stochastic schemes in geophysics. SIAM Journal on Numerical Analysis 42, 2257–2276.

Franzke, C.E., O’Kane, T.J., Berner, J., Williams, P.D., Lucarini, V., 2015. Stochastic climate theory and modeling. Wiley Interdisciplinary Reviews: Climate Change 6, 63–78.

Gent, P.R., McWilliams, J.C., 1990. Isopycnal mixing in ocean circulation models. Journal of Physical Oceanography 20, 150–155.

Gent, P.R., Willebrand, J., McDougall, T.J., McWilliams, J.C., 1995. Parameterising eddy-induced tracer transports in ocean circulation models. Journal of Physical Oceanography 25, 463–474.

Giannakis, D., 2019. Data-driven spectral decomposition and forecasting of ergodic dynamical system. Applied and Computational Harmonic Analysis 47, 338–396.

Griffies, S.M., 1998. The Gent-McWilliams skew flux. Journal of Physical Oceanography 28, 831–841.

Griffies, S.M., Hallberg, R.W., 2000. Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models. Monthly Weather Review 128, 2935–2946.

Grooms, I., Majda, A.J., 2014. Stochastic superparameterization in quasigeostrophic turbulence. Journal of Computational Physics 271, 78–98.
Grooms, I., Majda, A.J., Smith, K.S., 2015. Stochastic superparameterization in a quasigeostrophic model of the Antarctic Circumpolar Current. Ocean Modelling 85, 1-15.

Gugole, F., Franzke, C.E., 2019. Numerical development and evaluation of an energy conserving conceptual stochastic climate model. Mathematics of Climate and Weather Forecasting 5, 45-64.

Gugole, F., Franzke, C.E., 2020. Spatial covariance modeling for stochastic subgrid-scale parameterizations using dynamic mode decomposition. Journal of Advances in Modeling Earth Systems 12.

Haidvogel, D.B., McWilliams, J.C., Gent, P.R., 1992. Boundary current separation in a quasi-geostrophic, eddy-resolving ocean circulation model. Journal of Physical Oceanography 22, 882-902.

Hogg, A.M., Dewar, W.K., Killworth, P.D., Blundell, J.R., 2003. A quasi-geostrophic coupled model (Q-GCM). Monthly Weather Review 131, 2261-2278.

Hogg, A.M., Killworth, P.D., Blundell, J.R., 2004. Mechanisms of decadal variability of the wind-driven ocean circulation. Journal of Physical Oceanography 35.

Holm, D.D., 2015. Variational principles for stochastic fluid dynamics. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 471, 20140963.

Jansen, M.F., Adcroft, A., Khani, S., Kong, H., 2019. Toward an energetically consistent, resolution aware parameterization of ocean mesoscale eddies. Journal of Advances in Modeling Earth Systems 11, 2844-2860.

Jansen, M.F., Held, I.M., 2014. Parameterizing subgrid-scale eddy effects using energetically consistent backscatter. Ocean Modelling 80, 36-48.

Karspeck, A.R., Yeager, S., Danabasoglu, G., Hoar, T., Collins, N., Raeder, K., Anderson, J., Tribbia, J., 2013. An ensemble adjustment Kalman filter for the CCSM4 ocean component. Journal of Climate 26, 7392-7413.

Kjellsson, J., Zanna, L., 2017. The impact of horizontal resolution on energy transfers in global ocean models. Fluids 2, 45.

Kunita, H., 1997. Stochastic flows and stochastic differential equations. volume 24 of *Cambridge Studies in Advanced Mathematics*. Cambridge University Press.

Leibovich, S., 1980. On wave-current interaction theories of Langmuir circulations. Journal of Fluid Mechanics 99, 715-724.

Leith, C.E., 1971. Atmospheric predictability and two-dimensional turbulence. Journal of the Atmospheric Sciences 28, 145-161.

Li, L., 2021. Stochastic modeling and numerical simulation of ocean dynamics. Theories. Université Rennes 1. URL: https://hal.archives-ouvertes.fr/tel-03207741/file/PhD-thesis-Li2021.pdf.
Li, L., Mémin, E., Tissot, G., 2022. Stochastic parameterization with dynamic mode decomposition. preprint hal-03597550 URL: https://hal.archives-ouvertes.fr/hal-03597550.

Marshall, D.P., Maddison, J.R., Berloff, P., 2012. A framework for parameterizing eddy potential vorticity fluxes. Journal of Physical Oceanography 42, 539–548.

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

McWilliams, J.C., 2006. Fundamentals of Geophysical Fluid Dynamics. Cambridge University Press.

Mémin, E., 2014. Fluid flow dynamics under location uncertainty. Geophysical & Astrophysical Fluid Dynamics 108, 119–146.

Mezić, I., 2013. Analysis of fluid flows via spectral properties of the Koopman operator. Annual Review of Fluid Mechanics 45, 357–378.

Porta Mana, P., Zanna, L., 2014. Toward a stochastic parameterization of ocean mesoscale eddies. Ocean Modelling 79, 1–20.

Redi, M.H., 1982. Oceanic isopycnal mixing by coordinate rotation. Journal of Physical Oceanography 12, 1154–1158.

Resseguier, V., Li, L., Jouan, G., Derian, P., Mémin, E., Chapron, B., 2020. New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics. Archives of Computational Methods in Engineering, 1886–1784.

Resseguier, V., Mémin, E., Chapron, B., 2017a. Geophysical flows under location uncertainty, part I: Random transport and general models. Geophysical & Astrophysical Fluid Dynamics 111, 149–176.

Resseguier, V., Mémin, E., Chapron, B., 2017b. Geophysical flows under location uncertainty, part II: Quasi-geostrophic models and efficient ensemble spreading. Geophysical & Astrophysical Fluid Dynamics 111, 177–208.

Resseguier, V., Mémin, E., Chapron, B., 2017c. Geophysical flows under location uncertainty, part III: SQG and frontal dynamics under strong turbulence. Geophysical & Astrophysical Fluid Dynamics 111, 209–227.

Schmid, P., 2010. Dynamic mode decomposition of numerical and experimental data. Journal of Fluid Mechanics 656, 5–28.

Sirovich, L., 1987. Turbulence and the dynamics of coherent structures, part I: Coherent structures. Quarterly of Applied Mathematics 45, 561–571.

Smagorinsky, J., 1963. General circulation experiments with the primitive equations. Monthly Weather Review 91, 99–164.

Vallis, G.K., 2017. Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation. 2 ed., Cambridge University Press.
Zanna, L., Porta 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.