#  Stochastic data-driven parameterization of unresolved mesoscale eddies

Long Li111Centre Inria Rennes – Bretagne Atlantique, France Bruno
Deremble222Institut des Géosciences de l’Environnement, Grenoble, France Noé
Lahaye††footnotemark:  Etienne Mémin††footnotemark:

###### 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).

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 ($\boldsymbol{X}_{t}$) under
LU is described by the following stochastic differential equation (SDE):

$\mathrm{d}\boldsymbol{X}_{t}(\boldsymbol{x})=\boldsymbol{v}\big{(}\boldsymbol{X}_{t}(\boldsymbol{x}),t\big{)}\,\mathrm{d}t+\boldsymbol{\sigma}\big{(}\boldsymbol{X}_{t}(\boldsymbol{x}),t\big{)}\,\mathrm{d}\boldsymbol{B}_{t},\
\quad\boldsymbol{X}_{0}(\boldsymbol{x})=\boldsymbol{x}\in\mathcal{D},$ (2.1)

where $\boldsymbol{v}$ denotes the time-smooth resolved velocity that is both
spatially and temporally correlated,
$\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{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, $\\{\boldsymbol{B}_{t}\\}_{\scriptscriptstyle 0\leq t\leq T}$
is an $\boldsymbol{\mathrm{I}}_{d}$-cylindrical Brownian motion (Da Prato and
Zabczyk, 2014) on a filtered probability space
$(\Omega,\mathcal{F},\\{\mathcal{F}_{t}\\}_{\scriptscriptstyle 0\leq t\leq
T},\mathbb{P})$333$\Omega$ is a set of samples, $\mathcal{F}$ is a
$\sigma$-field (a collection of subsets of $\Omega$),
$\\{\mathcal{F}_{t}\\}_{\scriptscriptstyle 0\leq t\leq 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. 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,\ldots,d}$.

The spatial structure of the unresolved flow component is modeled by the
correlation operator, $\boldsymbol{\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]$,
$\boldsymbol{\sigma}(\boldsymbol{\cdot},t)[\bullet]$ is assumed to be a
Hilbert-Schmidt integral operator on $H$ with a bounded matrix kernel
$\breve{\boldsymbol{\sigma}}=(\breve{\sigma}_{ij})_{i,j=1,\ldots,d}$ such that

$\boldsymbol{\sigma}(\boldsymbol{x},t)\,\boldsymbol{f}=\int_{\mathcal{D}}\breve{\boldsymbol{\sigma}}(\boldsymbol{x},\boldsymbol{y},t)\boldsymbol{f}(\boldsymbol{y})\,\mathrm{d}\boldsymbol{y},\
\quad\boldsymbol{f}\in H,\ \quad\boldsymbol{x}\in\mathcal{D}.$ (2.2)

The composition of $\boldsymbol{\sigma}(\boldsymbol{\cdot},t)[\bullet]$ and
its adjoint $\boldsymbol{\sigma}^{*}(\boldsymbol{\cdot},t)[\bullet]$ is a
trace class operator on $H$ and admits eigenfunctions
$\boldsymbol{\xi}_{n}(\boldsymbol{\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:

$\boldsymbol{\sigma}(\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{B}_{t}=\sum_{n\in\mathbb{N}}\lambda_{n}^{1/2}(t)\boldsymbol{\xi}_{n}(\boldsymbol{x},t)\,\mathrm{d}\beta_{n}(t),$
(2.3)

where $\beta_{n}$ are independent (one-dimensional) standard Brownian motions.

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

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

Moreover, the point-wise ($\boldsymbol{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
$\boldsymbol{\sigma}\boldsymbol{\sigma}^{*}$, and denoted by
$\boldsymbol{a}(\boldsymbol{x},t)$, namely

$\boldsymbol{a}(\boldsymbol{x},t):=\int_{\mathcal{D}}\breve{\boldsymbol{\sigma}}(\boldsymbol{x},\boldsymbol{y},t)\breve{\boldsymbol{\sigma}}^{\scriptscriptstyle
T}(\boldsymbol{y},\boldsymbol{x},t)\,\mathrm{d}\boldsymbol{y}=\sum_{n\in\mathbb{N}}\lambda_{n}(t)\big{(}\boldsymbol{\xi}_{n}\boldsymbol{\xi}_{n}^{\scriptscriptstyle
T}\big{)}(\boldsymbol{x},t).$ (2.5)

As shown in Bauer et al. (2020a), the process
$\int_{0}^{t}\boldsymbol{a}(\boldsymbol{x},s)\,\mathrm{d}s$ actually
corresponds to the quadratic variation (Da Prato and Zabczyk, 2014) of
$\int_{0}^{t}\boldsymbol{\sigma}(\boldsymbol{x},s)\,\mathrm{d}\boldsymbol{B}_{s}$,
and it is a continuous random process of finite variation (hence time
differentiable). In the particular case of a nonrandom $\boldsymbol{\sigma}$,
$\boldsymbol{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 $\boldsymbol{a}$ represents the (possibly random) friction
coefficients of the unresolved fluid motions and the eigenvalues $\lambda_{n}$
have the unit of m2/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(\boldsymbol{X}_{t+\delta t},t+\delta t)=\Theta(\boldsymbol{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):

$\mathrm{\mathbb{D}}_{t}\Theta=\mathrm{d}_{t}\Theta+(\boldsymbol{v}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})\boldsymbol{\cdot}\boldsymbol{\nabla}\Theta-\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\nabla}\Theta)\,\mathrm{d}t=0,$
(2.6)

in which $\mathrm{\mathbb{D}}_{t}$ is introduced as a stochastic transport
operator and
$\mathrm{d}_{t}\Theta(\boldsymbol{x}):=\Theta(\boldsymbol{x},t+\delta
t)-\Theta(\boldsymbol{x},t)$ stands for the (forward) time-increment of
$\Theta$ at a fixed point $\boldsymbol{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
$\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}\boldsymbol{\cdot}\boldsymbol{\nabla}\Theta$
continuously backscatters random energy to the system through the quadratic
variation
$\nicefrac{{1}}{{2}}\,(\boldsymbol{\nabla}\Theta)\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\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 $\boldsymbol{v}^{\star}$ in (2.6) which is defined as

$\boldsymbol{v}^{\star}=\boldsymbol{v}-\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{a}+\boldsymbol{\sigma}^{*}(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}).$
(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
$\nicefrac{{1}}{{2}}\,\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{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
$\mathrm{\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\geq 1)$ of the tracer is materially and integrally
invariant, namely

$\mathrm{\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}\boldsymbol{x}=0.$
(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:

$\displaystyle\underline{\text{Momentum equations}}$
$\displaystyle\mathrm{\mathbb{D}}_{t}\boldsymbol{v}+f\boldsymbol{k}\times(\boldsymbol{u}\,\mathrm{d}t+\boldsymbol{\sigma}_{h}\mathrm{d}\boldsymbol{B}_{t})=(b\,\mathrm{d}t+\mathrm{d}b_{t}^{\sigma})\boldsymbol{k}-\boldsymbol{\nabla}(p\,\mathrm{d}t+\mathrm{d}p_{t}^{\sigma}),$
(2.9a) $\displaystyle\underline{\text{Continuity equations}}$
$\displaystyle\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}^{\star}=0,\
\quad\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}=0,$
(2.9b) $\displaystyle\underline{\text{Thermodynamic equation}}$
$\displaystyle\mathrm{\mathbb{D}}_{t}b+N^{2}(w^{\star}\,\mathrm{d}t+\sigma_{z}\mathrm{d}B_{t})=\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}_{\boldsymbol{\cdot}z}N^{2})\,\mathrm{d}t,$
(2.9c)

where $\boldsymbol{v}=(\boldsymbol{u},w)^{\scriptscriptstyle T}$,
$\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}=(\boldsymbol{\sigma}_{h}\mathrm{d}\boldsymbol{B}_{t},\sigma_{z}\mathrm{d}B_{t})^{\scriptscriptstyle
T}$, $\boldsymbol{k}=(0,0,1)^{\scriptscriptstyle 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}\overline{\rho}(z)/\rho_{0}$ is the Brunt-Väisälä
frequency with $g$ the gravity value and $\overline{\rho}$ the stationary
density in equilibrium that only depend on height,
$b(\boldsymbol{x},t)=-g\rho^{\prime}(\boldsymbol{x},t)/\rho_{0}$ is the
resolved buoyancy associated with the density anomaly $\rho^{\prime}$,
$\mathrm{d}b_{t}^{\sigma}/\mathrm{d}t$ is a zero-mean noise that models the
unresolved buoyancy fluctuations and $\boldsymbol{a}_{\boldsymbol{\cdot}z}$
denotes the $z$-column vector of tensor $\boldsymbol{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 $\boldsymbol{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}\boldsymbol{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=\mathcal{O}(1)$ or
$\epsilon\,\mbox{R}_{\mbox{o}}=\mathcal{O}(\mbox{R}_{\mbox{o}})$ with a small
Rossby number $\mbox{R}_{\mbox{o}}=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}\mathrm{d}B_{t}/\boldsymbol{\sigma}_{h}\mathrm{d}\boldsymbol{B}_{t}\sim(H/L)(\mbox{R}_{\mbox{o}}/\mbox{B}_{\mbox{u}})$
with the height scale $H\ll L$ and the Burger number
$\mbox{B}_{\mbox{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
$\mathrm{d}p_{t}^{\sigma}\sim\Phi=\epsilon^{1/2}f_{0}L^{2}$ and
$\mathrm{d}b_{t}^{\sigma}\sim\Phi/H$. All the resolved variables
($\boldsymbol{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 ($\boldsymbol{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:

$\displaystyle\underline{\text{Momentum equations}}$
$\displaystyle\mathrm{\mathbb{D}}_{t}^{h}\boldsymbol{u}+f_{0}\boldsymbol{k}\times\boldsymbol{u}_{a}\,\mathrm{d}t+\beta
y\boldsymbol{k}\times(\boldsymbol{u}\,\mathrm{d}t+\boldsymbol{\sigma}_{h}\mathrm{d}\boldsymbol{B}_{t})=-\boldsymbol{\nabla}_{h}p_{a}\,\mathrm{d}t,$
(2.10a) $\displaystyle\underline{\text{Buoyancy equation}}$
$\displaystyle\mathrm{\mathbb{D}}_{t}^{h}b+N^{2}w_{a}\,\mathrm{d}t=0,$ (2.10b)
$\displaystyle\underline{\text{Continuity equation}}$
$\displaystyle\boldsymbol{\nabla}_{h}\boldsymbol{\cdot}\boldsymbol{u}_{a}+\partial_{z}w_{a}=0,$
(2.10c) where $\boldsymbol{u}$ and $p$ denote the geostrophic components of
the resolved velocity and pressure,
$(\boldsymbol{u}_{a},w_{a})^{\scriptscriptstyle T}$ and $p_{a}$ denote the
ageostrophic components, and $\mathrm{\mathbb{D}}_{t}^{h}$ stands for the
horizontal component of the stochastic transport operator, which is defined as
$\mathrm{\mathbb{D}}_{t}^{h}\Theta=\mathrm{d}_{t}\Theta+(\boldsymbol{u}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}_{h}\mathrm{d}\boldsymbol{B}_{t})\boldsymbol{\cdot}\boldsymbol{\nabla}_{h}\Theta-\frac{1}{2}\,\boldsymbol{\nabla}_{h}\boldsymbol{\cdot}(\boldsymbol{a}_{h}\boldsymbol{\nabla}_{h}\Theta)\,\mathrm{d}t$
(2.10d) with
$\boldsymbol{u}^{\star}:=\boldsymbol{u}-\nicefrac{{1}}{{2}}\,\boldsymbol{\nabla}_{h}\boldsymbol{\cdot}\boldsymbol{a}_{h}$
and $\boldsymbol{\nabla}_{h}=(\partial_{x},\partial_{y})^{\scriptscriptstyle
T}$. The geostrophic component and the noise satisfy the following
equilibrium: $\displaystyle\underline{\text{Geostrophic balances and
martingale pressure equilibrium}}$
$\displaystyle\boldsymbol{u}=\frac{1}{f_{0}}\boldsymbol{\nabla}^{\perp}_{h}p,\
\quad\boldsymbol{\sigma}_{h}\mathrm{d}\boldsymbol{B}_{t}=\frac{1}{f_{0}}\boldsymbol{\nabla}^{\perp}_{h}\mathrm{d}p_{t}^{\sigma},\
\quad\boldsymbol{\nabla}_{h}\boldsymbol{\cdot}(\boldsymbol{u}^{\star}-\boldsymbol{u})=0,$
(2.10e) $\displaystyle\underline{\text{Hydrostatic balances}}$ $\displaystyle
b=\partial_{z}p,\
\quad\mathrm{d}b_{t}^{\sigma}=\partial_{z}\mathrm{d}p_{t}^{\sigma},$ (2.10f)

where
$\boldsymbol{\nabla}^{\perp}_{h}=(-\partial_{y},\partial_{x})^{\scriptscriptstyle
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}\boldsymbol{u}\perp\boldsymbol{\nabla}b$ and
$\partial_{z}\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}\perp\boldsymbol{\nabla}\mathrm{d}b_{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

$\displaystyle\underline{\text{Stochastic evolution of PV}}$
$\displaystyle\mathrm{\mathbb{D}}_{t}q=\sum_{i=1,2}-\mathrm{J}\big{(}(u^{\star}\,\mathrm{d}t+\sigma\mathrm{d}B_{t})^{i},u^{i}\big{)}+\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\big{(}\partial_{x_{i}}^{\scriptscriptstyle\perp}\boldsymbol{a}\boldsymbol{\nabla}u^{i}\big{)}\,\mathrm{d}t-\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\nabla}f)\,\mathrm{d}t$
$\displaystyle\>\phantom{\mathrm{\mathbb{D}}_{t}=}\>-\frac{f_{0}}{N^{2}}\partial_{z}(\boldsymbol{u}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})\boldsymbol{\cdot}\boldsymbol{\nabla}b+\frac{f_{0}}{2N^{2}}\boldsymbol{\nabla}\boldsymbol{\cdot}(\partial_{z}\boldsymbol{a}\boldsymbol{\nabla}b)\,\mathrm{d}t,$
(2.11a) $\displaystyle\underline{\text{From PV to streamfunction}}$
$\displaystyle
q=\nabla^{2}\psi+\partial_{z}\Big{(}\frac{f_{0}}{N^{2}}b\Big{)}+\beta y,\
\quad b=f_{0}\partial_{z}\psi$ (2.11b)
$\displaystyle\underline{\text{Incompressible constraints}}$
$\displaystyle\boldsymbol{u}=\boldsymbol{\nabla}^{\perp}\psi,\
\quad\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}=\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}^{\star}-\boldsymbol{u})=0,$
(2.11c)

where $q$ denotes the PV, $\psi:=p/f_{0}$ stands for the streamfunction of the
resolved flow,
$\mathrm{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, $\mathrm{J}(u^{i},u^{i})=0$, and to the
thermal wind balance, $\partial_{z}\boldsymbol{u}\perp\boldsymbol{\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

$\mathrm{d}_{t}\int_{\mathcal{D}}\frac{1}{2}\Big{(}|\boldsymbol{u}|^{2}+\frac{b^{2}}{N^{2}}\Big{)}\,\mathrm{d}\boldsymbol{x}=\Big{(}\int_{\mathcal{D}}w_{a}b\,\mathrm{d}\boldsymbol{x}-\int_{\mathcal{D}}w_{a}b\,\mathrm{d}\boldsymbol{x}\Big{)}\,\mathrm{d}t=0.$
(2.12)

The demonstration can be found in Appendix A. Conservation of path-wise
energy. 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 Appendix B. Conversion of ensemble energies. Based
on the decompositions
$\boldsymbol{u}=\mathbb{E}[\boldsymbol{u}]+(\boldsymbol{u}-\mathbb{E}[\boldsymbol{u}]):=\overline{\boldsymbol{u}}+\boldsymbol{u}^{\prime}$
and $b=\overline{b}+b^{\prime}$, 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
$\nicefrac{{1}}{{2}}\langle|\overline{\boldsymbol{u}}|^{2}\rangle$,
$\nicefrac{{1}}{{2}}\langle\overline{|\boldsymbol{u}^{\prime}|^{2}}\rangle$,
$\nicefrac{{1}}{{2}}\langle(\overline{b}/N)^{2}\rangle$ and
$\nicefrac{{1}}{{2}}\langle\overline{(b^{\prime}/N)^{2}}\rangle$, where
$\langle\boldsymbol{\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.

MKEEKEMPEEPE$-\,\big{\langle}\big{[}\boldsymbol{\nabla}\overline{\boldsymbol{u}},\overline{(\boldsymbol{u}^{\star})^{\prime}\boldsymbol{u}^{\prime}}\big{]}_{\scriptscriptstyle
F}\big{\rangle}+\nicefrac{{1}}{{2}}\big{\langle}\big{[}\boldsymbol{\nabla}\overline{\boldsymbol{u}},\
\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{\boldsymbol{u}}+\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\boldsymbol{u}^{\prime}}\big{]}_{\scriptscriptstyle
F}\big{\rangle}$$+\,\big{\langle}\big{[}\boldsymbol{\nabla}\overline{\boldsymbol{u}},\overline{(\boldsymbol{u}^{\star})^{\prime}\boldsymbol{u}^{\prime}}\big{]}_{\scriptscriptstyle
F}\big{\rangle}-\nicefrac{{1}}{{2}}\big{\langle}\big{[}\boldsymbol{\nabla}\overline{\boldsymbol{u}},\
\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{\boldsymbol{u}}+\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\boldsymbol{u}^{\prime}}\big{]}_{\scriptscriptstyle
F}\big{\rangle}$$-\,\big{\langle}(\boldsymbol{\nabla}\overline{b}/N^{2})\boldsymbol{\cdot}\overline{(\boldsymbol{u}^{\star})^{\prime}b^{\prime}}\big{\rangle}+\nicefrac{{1}}{{2}}\big{\langle}(\boldsymbol{\nabla}\overline{b}/N^{2})\boldsymbol{\cdot}(\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{b}+\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}b^{\prime}})\big{\rangle}$$+\,\big{\langle}(\boldsymbol{\nabla}\overline{b}/N^{2})\boldsymbol{\cdot}\overline{(\boldsymbol{u}^{\star})^{\prime}b^{\prime}}\big{\rangle}-\nicefrac{{1}}{{2}}\big{\langle}(\boldsymbol{\nabla}\overline{b}/N^{2})\boldsymbol{\cdot}(\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{b}+\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}b^{\prime}})\big{\rangle}$$-\,\big{\langle}\overline{b}\overline{w}_{a}\big{\rangle}$$+\,\big{\langle}\overline{b}\overline{w}_{a}\big{\rangle}$$-\,\big{\langle}\overline{b^{\prime}w_{a}^{\prime}}\big{\rangle}$$+\,\big{\langle}\overline{b^{\prime}w_{a}^{\prime}}\big{\rangle}$

Here, $[\boldsymbol{A},\boldsymbol{B}]_{\scriptscriptstyle
F}=\mathrm{Tr}(\boldsymbol{A}^{\scriptscriptstyle T}\boldsymbol{B})$ stands
for the Frobenius inner product of the tensors $\boldsymbol{A}$ and
$\boldsymbol{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 $\boldsymbol{\sigma}$ random, the
classical geostrophic Reynolds stress (Bachman, 2019) is modified to
$\overline{(\boldsymbol{u}^{\star})^{\prime}\boldsymbol{u}^{\prime}}$ 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
$\overline{\boldsymbol{a}}\,\mathrm{d}t=\mathrm{Var}(\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})$
is symmetric and non-negative definite, ensures that the integral
$\big{\langle}\big{[}\boldsymbol{\nabla}\overline{\boldsymbol{u}},\
\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{\boldsymbol{u}}\big{]}_{\scriptscriptstyle
F}\big{\rangle}$ (resp. $\big{\langle}\big{[}\boldsymbol{\nabla}\overline{b},\
\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{b}\big{]}_{\scriptscriptstyle
F}\big{\rangle}$) is always positive, hence provides a positive transfer from
MKE to EKE (resp. from MPE to EPE). When the correlation tensor
$\boldsymbol{\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 $\boldsymbol{\Gamma}_{t}$ be an $H$-valued $\mathcal{F}_{t}$-predictable
process satisfying
$\mathbb{E}_{\scriptscriptstyle\mathbb{P}}\big{[}\exp(\frac{1}{2}\int_{0}^{\scriptscriptstyle
T}\|\boldsymbol{\Gamma}_{t}\|_{\scriptscriptstyle
H}^{2}\,\mathrm{d}t)\big{]}<+\infty$, then the process
$\\{\widetilde{\boldsymbol{B}}_{t}:=\boldsymbol{B}_{t}+\int_{0}^{t}\boldsymbol{\Gamma}_{s}\,\mathrm{d}s\\}_{\scriptscriptstyle
0\leq t\leq T}$ is an $H$-valued cylindrical Wiener process on
$(\Omega,\mathcal{F},\\{\mathcal{F}_{t}\\}_{\scriptscriptstyle 0\leq t\leq
T},\widetilde{\mathbb{P}})$ by a change of the probability measure:

$\frac{\mathrm{d}\widetilde{\mathbb{P}}}{\mathrm{d}\mathbb{P}}=\exp\Big{(}-\int_{0}^{T}\left\langle{\boldsymbol{\Gamma}_{t},\mathrm{d}\boldsymbol{B}_{t}}\right\rangle_{\scriptscriptstyle
H}-\frac{1}{2}\int_{0}^{T}\|\boldsymbol{\Gamma}_{t}\|_{\scriptscriptstyle
H}^{2}\,\mathrm{d}t\Big{)},$ (2.13)

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

$\widetilde{\boldsymbol{u}^{\star}}:=\boldsymbol{u}^{\star}-\boldsymbol{\sigma}\boldsymbol{\Gamma}_{t}.$
(2.14)

Hereafter, the resulting stochastic system will be adopted with
$\boldsymbol{\sigma}\boldsymbol{\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},\boldsymbol{u}_{k},\boldsymbol{\sigma}_{k},\boldsymbol{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,\ldots,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}^{\prime}}{(H_{k}+H_{k+1})/2},$ (3.1)

where $g_{k}^{\prime}:=-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}}\Big{(}\frac{\psi_{k-1}-\psi_{k}}{g_{k-1}^{\prime}}-\frac{\psi_{k}-\psi_{k+1}}{g_{k}^{\prime}}\Big{)}=q_{k}-\beta
y.$ (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 $\boldsymbol{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

$\widetilde{\mathrm{\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)\,\mathrm{d}t,$
(3.2b)

where $\widetilde{\mathrm{\mathbb{D}}}_{t}^{k}$ denotes the $k$-th layer
stochastic transport operator under the probability measure
$\widetilde{\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 a conservative flux
form $\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}q)$ that corresponds
exactly to the 9-points Arakawa Jacobian scheme (Arakawa and Lamb, 1981) of
$\mathrm{J}(\psi,q)$, in which the zonal and meridional components of the
advection flux are defined as
$\nicefrac{{2}}{{3}}\,\big{(}\overline{\overline{u}^{x}}^{y}\overline{q}^{x}\big{)}+\nicefrac{{1}}{{3}}\,\big{(}\overline{\overline{u}^{x}\overline{q}^{d}}^{y}\big{)}$
and
$\nicefrac{{2}}{{3}}\,\big{(}\overline{\overline{v}^{y}}^{x}\overline{q}^{y}\big{)}+\nicefrac{{1}}{{3}}\,\big{(}\overline{\overline{v}^{x}\overline{q}^{d}}^{y}\big{)}$.
Here, $\overline{\boldsymbol{\cdot}}^{x},\overline{\boldsymbol{\cdot}}^{y}$
and $\overline{\boldsymbol{\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
$\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\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}=\boldsymbol{\nabla}^{\perp}\boldsymbol{\cdot}\boldsymbol{\tau}/f_{0}$
and $e_{n}=\nicefrac{{1}}{{2}}\,\delta_{ek}\nabla^{2}\psi_{n}$, where
$\boldsymbol{\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.

Figure 1: Illustration for the vertical discretization of a three-layered QG
system.

###  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 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
$\\{\boldsymbol{u}_{\scriptscriptstyle\text{HR}}(\boldsymbol{x},t_{i})\\}_{i=1,\ldots,N}$
be the finite set of velocity snapshots provided by a high-resolution (HR)
simulation. We first build the spatial local fluctuations
$\boldsymbol{u}_{f}(\boldsymbol{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_{\scriptscriptstyle\text{HR}}$, to obtain the streamfunction
fluctuations,
$\psi_{f}(\boldsymbol{x},t_{i})=\big{(}(\mathrm{I}-\mathrm{G})\star\psi_{\scriptscriptstyle\text{HR}}\big{)}(\boldsymbol{x},t_{i})$
(only for the coarse-grid points $\boldsymbol{x}$). Then, the geostrophic
velocity fluctuations can be derived by
$\boldsymbol{u}_{f}=\boldsymbol{\nabla}^{\perp}_{\scriptscriptstyle\text{LR}}\psi_{f}$.

We next center the data set by
$\boldsymbol{u}_{f}^{\prime}=\boldsymbol{u}_{f}-\overline{\boldsymbol{u}_{f}}^{t}$
(with $\overline{\bullet}^{t}$ the temporal mean) and perform the EOF
procedure (Sirovich, 1987) to get a set of orthogonal temporal modes
$\\{\alpha_{n}\\}_{n=1,\ldots,N}$ and orthonormal spatial modes
$\\{\boldsymbol{\phi}_{n}\\}_{n=1,\ldots,N}$ satisfying

$\boldsymbol{u}_{f}^{\prime}(\boldsymbol{x},t_{i})=\sum_{n=1}^{N}\alpha_{n}(t_{i})\boldsymbol{\phi}_{n}(\boldsymbol{x}),$
(3.3)

where
$\langle\boldsymbol{\phi}_{m},\boldsymbol{\phi}_{n}\rangle_{\scriptscriptstyle\mathcal{H}}=\delta_{m,n}$,
$\overline{\alpha_{m}\alpha_{n}}^{t}=\lambda_{n}\delta_{m,n}$,
$\mathcal{H}:=L^{2}(\mathcal{D})$ and $\mathcal{D}=\mathcal{A}\times[-H,0]$
denotes the 3D domain with the horizontal area $\mathcal{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

$\displaystyle\boldsymbol{\sigma}_{k}(\boldsymbol{x})\,\mathrm{d}\widetilde{\boldsymbol{B}}_{t}$
$\displaystyle=\tau^{1/2}\sum_{n=1}^{M}\lambda_{n}^{1/2}\boldsymbol{\phi}_{n,k}(\boldsymbol{x})\,\mathrm{d}\beta_{n}(t),$
(3.4a) $\displaystyle\boldsymbol{a}_{k}(\boldsymbol{x})$
$\displaystyle=\tau\sum_{n=1}^{M}\lambda_{n}\boldsymbol{\phi}_{n,k}(\boldsymbol{x})\boldsymbol{\phi}_{n,k}^{\scriptscriptstyle
T}(\boldsymbol{x}).$ (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_{\scriptscriptstyle\text{LR}}$.

In this work, we simply propose to represent the Girsanov drift by the
temporal mean of the unresolved eddies, that is

$\boldsymbol{\sigma}_{k}(\boldsymbol{x})\boldsymbol{\Gamma}_{t}=\sum_{n=1}^{N}\langle\overline{\boldsymbol{u}_{f}}^{t},\boldsymbol{\phi}_{n}\rangle_{\scriptscriptstyle\mathcal{H}}\boldsymbol{\phi}_{n,k}(\boldsymbol{x}).$
(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,
$\widetilde{\boldsymbol{\sigma}}\mathrm{d}\widetilde{\boldsymbol{B}}_{t}\perp\boldsymbol{\nabla}\partial_{z}\widetilde{b}$
with $\widetilde{b}:=f_{0}b/N^{2}$. This projection procedure reads

$\widetilde{\boldsymbol{\sigma}}_{k}(\boldsymbol{x},t)\,\mathrm{d}\widetilde{\boldsymbol{B}}_{t}=A^{1/2}\,\boldsymbol{\mathrm{P}}_{k}(\boldsymbol{x},t)\
\boldsymbol{\sigma}_{k}(\boldsymbol{x})\,\mathrm{d}\widetilde{\boldsymbol{B}}_{t},\
\quad\boldsymbol{\mathrm{P}}_{k}:=\boldsymbol{\mathrm{I}}_{2}-\frac{\boldsymbol{\nabla}(\partial_{z}\widetilde{b})_{k}\big{(}\boldsymbol{\nabla}(\partial_{z}\widetilde{b})_{k}\big{)}^{\scriptscriptstyle
T}}{|\boldsymbol{\nabla}(\partial_{z}\widetilde{b})_{k}|^{2}},$ (3.5)

where $\boldsymbol{\sigma}_{k}\mathrm{d}\widetilde{\boldsymbol{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 $\boldsymbol{a}$ defined through (3.4b) is updated as
$\widetilde{\boldsymbol{a}}_{k}=A\,\boldsymbol{\mathrm{P}}_{k}\,\boldsymbol{a}_{k}\,\boldsymbol{\mathrm{P}}_{k}^{\scriptscriptstyle
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 $\boldsymbol{\tau}=(-\tau_{0}\cos(2\pi
y)/L_{y},0)^{\scriptscriptstyle T}$. A mixed horizontal boundary condition is
used for the $k$-th layer streamfunction:
$\psi_{k}|_{\partial\mathcal{A}}=f_{k}(t)$ and
$\partial_{n}^{2}\psi_{k}|_{\partial\mathcal{A}}=-(\alpha_{\text{bc}}/\Delta
x)\partial_{n}\psi_{k}|_{\partial\mathcal{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_{\text{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.

Parameters | Value | Description
---|---|---
$X\times Y$ |  $(3840\times 4800)$ km | Domain size
$H_{k}$ |  $(350,750,2900)$ m | Mean layer thickness
$\rho$ |  $1000$ kg m3 | Density
$g_{k}^{\prime}$ |  $(0.025,0.0125)$ m s-2 | Reduced gravity
$\delta_{\text{ek}}$ |  $2$ m | Bottom Ekman layer thickness
$\tau_{0}$ |  $2\times 10^{-5}$ m2 s-2 | Wind stress magnitude
$\alpha_{\text{bc}}$ | $0.2$ | Mixed boundary condition coefficient
$f_{0}$ |  $9.375\times 10^{-5}$ s-1 | Mean Coriolis parameter
$\beta$ |  $1.754\times 10^{-11}$ (m s)-1 | Coriolis parameter gradient
$L_{d}$ |  $(39,22)$ km | Baroclinic Rossby radii

Table 1: Common parameters for all the models.

$\Delta x$ (km) |  $\Delta t$ (s) |  $A_{4}$ (m4 s-1)
---|---|---
5 | 600 | $2\times 10^{9}$
40 | 1200 | $5\times 10^{11}$
80 | 1440 | $5\times 10^{12}$
120 | 1800 | $1\times 10^{13}$

Table 2: Values of grid varying parameters.

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}}\boldsymbol{\phi}_{1}$) and the noise
variance ($\boldsymbol{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 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.

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.

###  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.

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}$ m3/s.

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 is defined as

$\text{GRE}=\frac{1}{|\mathcal{D}|}\int_{\mathcal{D}}\frac{1}{2}\left(\frac{\big{(}\overline{\psi^{\scriptscriptstyle\text{M}}}^{t}-\overline{\psi^{\scriptscriptstyle\text{R}}}^{t}\big{)}^{2}}{\sigma^{2}[\psi^{\scriptscriptstyle\text{M}}]}+\frac{\sigma^{2}[\psi^{\scriptscriptstyle\text{R}}]}{\sigma^{2}[\psi^{\scriptscriptstyle\text{M}}]}-1-\log\Big{(}\frac{\sigma^{2}[\psi^{\scriptscriptstyle\text{R}}]}{\sigma^{2}[\psi^{\scriptscriptstyle\text{M}}]}\Big{)}\right)\,\mathrm{d}\boldsymbol{x},$
(4.1)

where $\sigma$ denotes the std, $\psi^{\scriptscriptstyle\text{M}}$ (resp.
$\psi^{\scriptscriptstyle\text{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.

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.

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.

Figure 6: Temporal mean of vertically integrated KE spectra for different
models (by colors) at different resolutions (by sub-figures).

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 $\boldsymbol{u}$ and $b$ is simply build 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.

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/m3, EKE = 0.62 J/m3, MPE
= 4.01 J/m3 and EPE = 0.34 J/m3.

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.

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).

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}=-\nicefrac{{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

$\widetilde{\mathrm{\mathbb{D}}}_{t}^{m}T_{m}+T_{m}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{m}\,\mathrm{d}t=\left((K_{2}\nabla^{2}-K_{4}\nabla^{4})T_{m}+\frac{T_{m}+T_{1}}{2H_{m}}e_{0}\right)\,\mathrm{d}t,$
(4.2)

where
$\boldsymbol{u}_{m}=\boldsymbol{\nabla}^{\perp}\psi_{1}-\boldsymbol{\tau}^{\perp}/(f_{0}H_{m})$
is the divergent mixed layer velocity,
$\widetilde{\mathrm{\mathbb{D}}}_{t}^{m}$ is the mixed layer stochastic
transport operator under probability measure $\widetilde{\mathbb{P}}$
(associated to $\boldsymbol{u}_{m}$), $K_{2}=100$ m2/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
Leapfroq 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.

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.

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.

Figure 10: Comparison of time-mean contour of barotropic (top row) and first
baroclinic (bottom row) streamfunctions for different models (by columns).

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 11: Temporal mean of vertically integrated KE spectra for different
models (by colors) at different resolutions (by sub-figures).

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.

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/m3.

Figure 13: Comparison of the vertically-integrated EKE density (top row) and
EPE density (bottom row) for different models (by columns).

##  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 $\sim$ 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 hydrostasy.

Another possible research axis focuses on new parameterization methods for the
unresolved flow 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 $\mathrm{d}_{t}X=F\,\mathrm{d}t+G\,\mathrm{d}B_{t}$ and
$\mathrm{d}_{t}Y=F^{\prime}\,\mathrm{d}t+G^{\prime}\,\mathrm{d}B_{t}$ be two
S(P)DEs driven by a Brownian motion $B$ satisfying
$\mathbb{E}[G\mathrm{d}B_{t}]=\mathbb{E}[G^{\prime}\mathrm{d}B_{t}]=0$, then

$\displaystyle\mathrm{d}_{t}(XY)$
$\displaystyle=X\mathrm{d}_{t}Y+Y\mathrm{d}_{t}X+\mathrm{d}\langle
X,Y\rangle_{t}$
$\displaystyle=X\mathrm{d}_{t}Y+Y\mathrm{d}_{t}X+GG^{\prime}\,\mathrm{d}t,$
(5.1)

where $\langle
X,Y\rangle_{t}\stackrel{{\scriptstyle\mathbb{P}}}{{=}}\lim_{\scriptscriptstyle
n\rightarrow+\infty}\sum_{i=1}^{p_{n}}\big{(}X_{t_{i}}^{n}-X_{t_{i-1}}^{n}\big{)}\big{(}Y_{t_{i}}^{n}-Y_{t_{i-1}}^{n}\big{)}$
stands for the joint quadratic variation between the processes $X$ and $Y$
with $0=t_{0}^{n}<t_{1}^{n}<\cdots<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
$\boldsymbol{u}=(u,v)^{\scriptscriptstyle T}$ and equation (2.10b) for the
resolved buoyancy $b$,

$\displaystyle\mathrm{d}_{t}u$
$\displaystyle=-(\boldsymbol{u}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})\boldsymbol{\cdot}\boldsymbol{\nabla}u+\Big{(}\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\nabla}u)+f_{0}v_{a}+\beta
yv-\partial_{x}p_{a}\Big{)}\,\mathrm{d}t$ (5.2a)
$\displaystyle\mathrm{d}_{t}v$
$\displaystyle=-(\boldsymbol{u}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})\boldsymbol{\cdot}\boldsymbol{\nabla}v+\Big{(}\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\nabla}v)-f_{0}u_{a}-\beta
yu-\partial_{y}p_{a}\Big{)}\,\mathrm{d}t$ (5.2b)
$\displaystyle\mathrm{d}_{t}b$
$\displaystyle=-(\boldsymbol{u}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})\boldsymbol{\cdot}\boldsymbol{\nabla}b+\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\nabla}b)\,\mathrm{d}t-N^{2}w_{a}\,\mathrm{d}t,$
(5.2c)

together with the geostrophic equilibrium (2.10e), the hydrostatic balances
(2.10f) and the continuity equation (2.10c). Note that
$\boldsymbol{a}\,\mathrm{d}t$ corresponds to the quadratic variation of the
noise $\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}$, which is
$\boldsymbol{a}$ defined in (2.5). Applying equation (Appendix A. Conservation
of path-wise energy) 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

$\displaystyle\int_{\mathcal{D}}\mathrm{d}_{t}\Big{(}\frac{b^{2}}{2N^{2}}\Big{)}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=\int_{\mathcal{D}}\frac{1}{N^{2}}\Big{(}b\mathrm{d}_{t}b+\frac{1}{2}\mathrm{d}\langle
b,b\rangle_{t}\Big{)}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=\int_{\mathcal{D}}\frac{1}{N^{2}}\Big{(}-b(\boldsymbol{u}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})\boldsymbol{\cdot}\boldsymbol{\nabla}b+\frac{1}{2}b\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}\boldsymbol{\nabla}b)\,\mathrm{d}t$
$\displaystyle\>\phantom{=\int_{\mathcal{D}}\frac{1}{N^{2}}\Big{(}}\>-N^{2}w_{a}b\,\mathrm{d}t+\frac{1}{2}(\boldsymbol{\nabla}b)\boldsymbol{\cdot}\boldsymbol{a}\boldsymbol{\nabla}b\,\mathrm{d}t\Big{)}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=-\int_{\mathcal{D}}\boldsymbol{\nabla}\boldsymbol{\cdot}\Big{(}(\boldsymbol{u}^{\star}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t})\frac{b^{2}}{2N^{2}}\Big{)}\,\mathrm{d}\boldsymbol{x}-\int_{\mathcal{D}}w_{a}b\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle=-\int_{\mathcal{D}}w_{a}b\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$
(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}\boldsymbol{\cdot}\boldsymbol{n}\
|_{\partial\mathcal{A}}=\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}\boldsymbol{\cdot}\boldsymbol{n}\
|_{\partial\mathcal{A}}=0$ (which leads to
$\boldsymbol{u}^{\star}\boldsymbol{\cdot}\boldsymbol{n}\
|_{\partial\mathcal{A}}=0$ from equations (2.7) and (2.5)). Similarly to the
previous calculation, one can show that

$\displaystyle\int_{\mathcal{D}}\mathrm{d}_{t}\Big{(}\frac{u^{2}}{2}\Big{)}\,\mathrm{d}\boldsymbol{x}=\int_{\mathcal{D}}u(f_{0}v_{a}+\beta
yv-\partial_{x}p_{a})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$ (5.4a)
$\displaystyle\int_{\mathcal{D}}\mathrm{d}_{t}\Big{(}\frac{v^{2}}{2}\Big{)}\,\mathrm{d}\boldsymbol{x}=-\int_{\mathcal{D}}v(f_{0}u_{a}+\beta
yu+\partial_{y}p_{a})\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t.$ (5.4b)

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

$\displaystyle\int_{\mathcal{D}}\mathrm{d}_{t}\Big{(}\frac{|\boldsymbol{u}|^{2}}{2}\Big{)}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=-\int_{\mathcal{D}}\big{(}(\boldsymbol{\nabla}^{\perp}p)\boldsymbol{\cdot}\boldsymbol{u}_{a}^{\perp}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}p_{a})\big{)}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle=-\int_{\mathcal{D}}\big{(}\boldsymbol{\nabla}^{\perp}\boldsymbol{\cdot}(p\boldsymbol{u}_{a}^{\perp})-p\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_{a}\big{)}\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t$
$\displaystyle=-\int_{\mathcal{A}}\int_{-H}^{0}p\partial_{z}w_{a}\,\mathrm{d}z\,\mathrm{d}\mathcal{A}\,\mathrm{d}t$
$\displaystyle=-\int_{\mathcal{A}}\Big{(}\big{[}pw_{a}\big{]}_{\scriptscriptstyle-H}^{0}-\int_{-H}^{0}w_{a}\partial_{z}p\,\mathrm{d}z\Big{)}\,\mathrm{d}\mathcal{A}\,\mathrm{d}t$
$\displaystyle=\int_{\mathcal{D}}w_{a}b\,\mathrm{d}\boldsymbol{x}\,\mathrm{d}t,$
(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 $\boldsymbol{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 (Appendix A. Conservation of path-wise energy) and (Appendix
A. Conservation of path-wise energy), one show that the stochastic system
(5.2) conserves the path-wise total energy of the resolved geostrophic flow
component, namely

$\mathrm{d}_{t}\int_{\mathcal{D}}\frac{1}{2}\Big{(}|\boldsymbol{u}|^{2}+\frac{b^{2}}{N^{2}}\Big{)}\,\mathrm{d}\boldsymbol{x}=\int_{\mathcal{D}}\mathrm{d}_{t}\Big{(}\frac{|\boldsymbol{u}|^{2}}{2}\Big{)}\,\mathrm{d}\boldsymbol{x}+\int_{\mathcal{D}}\mathrm{d}_{t}\Big{(}\frac{b^{2}}{2N^{2}}\Big{)}\,\mathrm{d}\boldsymbol{x}=0.$
(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 $\boldsymbol{a}$) of (5.2)
into a mean component $\overline{\Theta}:=\mathbb{E}[\Theta]$ and an eddy
component $\Theta^{\prime}:=\Theta-\mathbb{E}[\Theta]$. The other time-smooth
components $F\,\mathrm{d}t$ apart from the geostrophic transport operator are
decomposed in the same way. For instance,
$F=-N^{2}(\overline{w_{a}}+w_{a}^{\prime})$ in the buoyancy equation. However,
the noise is only an eddy component since
$\mathbb{E}[\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}]=0$. Substituting
such decomposition into equations (5.2) and taking subsequently the ensemble-
mean of the expanded equation, we obtain

$\displaystyle\partial_{t}\overline{\Theta}=\overline{F}-\overline{\boldsymbol{u}^{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla}\overline{\Theta}+\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\big{(}\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{\Theta}\big{)}-\overline{(\boldsymbol{u}^{\star})^{\prime}\boldsymbol{\cdot}\boldsymbol{\nabla}\Theta^{\prime}}+\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}(\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\Theta^{\prime}}),$
(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, on can easily show the time evolution of the global energy for
$\overline{\Theta}$:
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{D}}\frac{\overline{\Theta}^{2}}{2}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=\int_{\mathcal{D}}\overline{F}\,\overline{\Theta}\,\mathrm{d}\boldsymbol{x}+\int_{\mathcal{D}}\overline{(\boldsymbol{u}^{\star})^{\prime}\Theta^{\prime}}\boldsymbol{\cdot}\boldsymbol{\nabla}\overline{\Theta}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle-\frac{1}{2}\int_{\mathcal{D}}\big{(}\boldsymbol{\nabla}\overline{\Theta}\big{)}\boldsymbol{\cdot}\big{(}\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{\Theta}+\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\Theta^{\prime}}\big{)}\,\mathrm{d}\boldsymbol{x}.$
(5.7b)

Applying this result for $\overline{u}$ and $\overline{v}$ (resp. for
$\overline{b}$) together with (Appendix A. Conservation of path-wise energy)
(resp. with (Appendix A. Conservation of path-wise energy)), 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^{\prime}$:

$\displaystyle\mathrm{d}_{t}\Theta^{\prime}=$
$\displaystyle-\Big{(}\overline{\boldsymbol{u}^{\star}}\boldsymbol{\cdot}\boldsymbol{\nabla}\Theta^{\prime}-\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\big{(}\overline{\boldsymbol{a}}\boldsymbol{\nabla}\Theta^{\prime}\big{)}\Big{)}\,\mathrm{d}t$
$\displaystyle-\big{(}(\boldsymbol{u}^{\star})^{\prime}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}\big{)}\boldsymbol{\cdot}\boldsymbol{\nabla}\overline{\Theta}-\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\big{(}\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\overline{\Theta}\big{)}\,\mathrm{d}t+F^{\prime}\,\mathrm{d}t$
$\displaystyle-\big{(}(\boldsymbol{u}^{\star})^{\prime}\,\mathrm{d}t+\boldsymbol{\sigma}\mathrm{d}\boldsymbol{B}_{t}\big{)}\boldsymbol{\cdot}\boldsymbol{\nabla}\Theta^{\prime}-\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}\big{(}\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\Theta^{\prime}\big{)}\,\mathrm{d}t$
$\displaystyle+\Big{(}\overline{(\boldsymbol{u}^{\star})^{\prime}\boldsymbol{\cdot}\boldsymbol{\nabla}\Theta^{\prime}-\frac{1}{2}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\Theta^{\prime})}\Big{)}\,\mathrm{d}t,$
(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^{\prime}$, the Itô’s integration-by-part formula (Appendix A.
Conservation of path-wise energy) is required again. In this case, one can
first show that the mean of the joint quadratic variation process for
$\Theta^{\prime}$ reads
$\displaystyle\overline{\mathrm{d}\langle\Theta^{\prime},\Theta^{\prime}\rangle_{t}}$
$\displaystyle=\overline{(\boldsymbol{\nabla}\Theta^{\prime})\boldsymbol{\cdot}(\overline{\boldsymbol{a}}\boldsymbol{\nabla}\Theta^{\prime})}+\overline{(\boldsymbol{\nabla}\Theta^{\prime})\boldsymbol{\cdot}(\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\Theta^{\prime})}+\overline{(\boldsymbol{\nabla}\Theta^{\prime})\boldsymbol{\cdot}(\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\overline{\Theta})}$
$\displaystyle+(\boldsymbol{\nabla}\overline{\Theta})\boldsymbol{\cdot}(\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\Theta^{\prime}})+(\boldsymbol{\nabla}\overline{\Theta})\boldsymbol{\cdot}(\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{\Theta}).$
(5.8b) Subsequently, we can deduce the evolution of the global variance for
$\Theta$, namely
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\int_{\mathcal{D}}\overline{\frac{(\Theta^{\prime})^{2}}{2}}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle=\int_{\mathcal{D}}\overline{F^{\prime}\Theta^{\prime}}\,\mathrm{d}\boldsymbol{x}-\int_{\mathcal{D}}\overline{(\boldsymbol{u}^{\star})^{\prime}\Theta^{\prime}}\boldsymbol{\cdot}\boldsymbol{\nabla}\overline{\Theta}\,\mathrm{d}\boldsymbol{x}$
$\displaystyle+\frac{1}{2}\int_{\mathcal{D}}\boldsymbol{\nabla}(\overline{\Theta})\boldsymbol{\cdot}(\overline{\boldsymbol{a}}\boldsymbol{\nabla}\overline{\Theta}+\overline{\boldsymbol{a}^{\prime}\boldsymbol{\nabla}\Theta^{\prime}})\,\mathrm{d}\boldsymbol{x}.$
(5.8c)

Applying similarly this result for $\overline{u}$ and $\overline{v}$ (resp.
for $\overline{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 et al. (2019) 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 and Lamb (1981) Arakawa, A., Lamb, V.R., 1981\.  A potential enstrophy and energy conserving scheme for the shallow water equations.  Monthly Weather Review 109, 18–36.
  * Arbic et al. (2013) Arbic, B.K., Polzin, K.L., Scott, R.B., Richman, J.G., Shriver, J.F., 2013.  On eddy viscosity, energy cascades, and the horizonal resolution of gridded stallite altimeter products.  Journal of Physical Oceanography 43, 283–300.
  * Bachman (2019) 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 et al. (2017) 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 et al. (2020a) 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 et al. (2020b) 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 (2005) Berloff, P., 2005.  Random-forcing model of the mesoscale oceanic eddies.  Journal of Fluid Mechanics 529, 71–95.
  * Berloff et al. (2007) 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 (2013) Brankart, J., 2013.  Impact of uncertainties in the horizontal density gradient upon low resolution global ocean modeling.  Ocean Modelling 66, 64–76.
  * Brecht et al. (2021) 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 et al. (2019) 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 et al. (2020) 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 and Leibovich (1976) Craik, A., Leibovich, S., 1976\.  A rational model for Langmuir circulations.  Journal of Fluid Mechanics 73, 401–426.
  * Da Prato and Zabczyk (2014) Da Prato, G., Zabczyk, J., 2014\.  Stochastic equations in infinite dimensions.  Encyclopedia of Mathematics and its Applications. 2 ed., Cambridge University Press.
  * Dufée et al. (2022) 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 and Témam (2005) Ewald, B., Témam, R., 2005\.  Numerical analysis of stochastic schemes in geophysics.  SIAM Journal on Numerical Analysis 42, 2257–2276.
  * Franzke et al. (2015) 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 and McWilliams (1990) Gent, P.R., McWilliams, J.C., 1990\.  Isopycnal mixing in ocean circulation models.  Journal of Physical Oceanography 20, 150–155.
  * Gent et al. (1995) 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 (2019) Giannakis, D., 2019.  Data-driven spectral decomposition and forecasting of ergodic dynamical system.  Applied and Computational Harmonic Analysis 47, 338–396.
  * Griffies (1998) Griffies, S.M., 1998.  The Gent-McWilliams skew flux.  Journal of Physical Oceanography 28, 831–841.
  * Griffies and Hallberg (2000) 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 and Majda (2014) Grooms, I., Majda, A.J., 2014\.  Stochastic superparameterization in quasigeostrophic turbulence.  Journal of Computational Physics 271, 78–98.
  * Grooms et al. (2015) 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 and Franzke (2019) 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 and Franzke (2020) 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 et al. (1992) Haidvogel, D.B., McWilliams, J.C., Gent, P.R., 1992.  Boundary current separation in a quasigeostrophic, eddy-resolving ocean circulation model.  Journal of Physical Oceanography 22, 882–902.
  * Hogg et al. (2003) 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 et al. (2004) 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 (2015) 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 et al. (2019) 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 and Held (2014) Jansen, M.F., Held, I.M., 2014\.  Parameterizing subgrid-scale eddy effects using energetically consistent backscatter.  Ocean Modelling 80, 36–48.
  * Karspeck et al. (2013) Karspeck, A.R., Yeager, S., Danabasoglu, G., Hoar, T., Collins, N., Raeder, K., Anderson, J., Tribbia, J., 2013\.  An ensemble ajustment Kalman filter for the CCSM4 ocean component.  Journal of Climate 26, 7392–7413.
  * Kjellsson and Zanna (2017) Kjellsson, J., Zanna, L., 2017\.  The impact of horizontal resolution on energy transfers in global ocean models.  Fluids 2, 45.
  * Kunita (1997) Kunita, H., 1997.  Stochastic flows and stochastic differential equations. volume 24 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press.
  * Leibovich (1980) Leibovich, S., 1980.  On wave-current interaction theories of Langmuir circulations.  Journal of Fluid Mechanics 99, 715–724.
  * Leith (1971) Leith, C.E., 1971.  Atmospheric predictability and two-dimensional turbulence.  Journal of the Atmospheric Sciences 28, 145–161.
  * Li (2021) Li, L., 2021.  Stochastic modeling and numerical simulation of ocean dynamics.  Theses. Université Rennes 1.  URL: https://hal.archives-ouvertes.fr/tel-03207741/file/PhD-thesis-Li2021.pdf.
  * Li et al. (2022) 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 et al. (2012) 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 (1977) 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 (2006) McWilliams, J.C., 2006.  Fundamentals of Geophysical Fluid Dynamics.  Cambridge University Press.
  * Mémin (2014) Mémin, E., 2014.  Fluid flow dynamics under location uncertainty.  Geophysical & Astrophysical Fluid Dynamics 108, 119–146.
  * Mezić (2013) Mezić, I., 2013.  Analysis of fluid flows via spectral properties of the Koopman operator.  Annual Review of Fluid Mechanics 45, 357–378.
  * Porta Mana and Zanna (2014) Porta Mana, P., Zanna, L., 2014\.  Toward a stochastic parameterization of ocean mesoscale eddies.  Ocean Modelling 79, 1–20.
  * Redi (1982) Redi, M.H., 1982.  Oceanic isopycnal mixing by coordinate rotation.  Journal of Physical Oceanography 12, 1154–1158.
  * Resseguier et al. (2020) 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 et al. (2017a) 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 et al. (2017b) 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 et al. (2017c) 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 (2010) Schmid, P., 2010.  Dynamic mode decomposition of numerical and experimental data.  Journal of Fluid Mechanics 656, 5–28.
  * Sirovich (1987) Sirovich, L., 1987.  Turbulence and the dynamics of coherent structures, part I: Coherent structures.  Quarterly of Applied Mathematics 45, 561–571.
  * Smagorinsky (1963) Smagorinsky, J., 1963.  General circulation experiments with the primitive equations.  Monthly Weather Review 91, 99–164.
  * Vallis (2017) Vallis, G.K., 2017.  Atmospheric and oceanic fluid dynamics: fundamentals and large-scale circulation.  2 ed., Cambridge University Press.
  * Zanna et al. (2017) 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.