Then,

$$
\begin{align*}
J &= d_t \left\langle \int_0^t \left( \frac{\partial}{\partial x_i} (d_t q \sigma_{s,ij}) \right) \cdot \int_0^t d\mathbf{B}_s^j \right\rangle \\
&= d_t \left\langle \int_0^t \left( \frac{\partial}{\partial x_i} \left( \left( \mathbf{Q}_\sigma \cdot d\mathbf{B}_s - \frac{\partial}{\partial x_k} (q\sigma_{s,lm} d\mathbf{B}_s^m) \right) \sigma_{s,ij} \right) \right) \cdot \int_0^t d\mathbf{B}_s^j \right\rangle \\
&= \frac{\partial}{\partial x_i} (\mathbf{Q}_\sigma^j \sigma_{t,ij}) \, dt - \frac{\partial}{\partial x_i} \left( \frac{\partial}{\partial x_l} (q\sigma_{t,lj}) \sigma_{t,ij} \right) \, dt \\
&= \frac{\partial}{\partial x_i} (\mathbf{Q}_\sigma^j \sigma_{t,ij}) \, dt - \frac{\partial}{\partial x_i} \left( \frac{\partial q}{\partial x_l} \sigma_{t,lj} \sigma_{t,ij} \right) \, dt - \frac{\partial}{\partial x_i} \left( q \frac{\partial \sigma_{t,lj}}{\partial x_l} \sigma_{t,ij} \right) \, dt \\
&= \nabla \cdot (\boldsymbol{\sigma}_t \mathbf{Q}_\sigma) \, dt - \nabla \cdot (\mathbf{a} \nabla q) \, dt - \nabla \cdot (q \boldsymbol{\sigma}_t (\nabla \cdot \boldsymbol{\sigma}_t)) \, dt. \tag{64}
\end{align*}
$$

In addition, we make the hypothesis that $d\mathbf{Q}_\sigma$ is time-differentiable in the Lagrangian frame, such that we have

$$
\begin{align}
d \int_{V(t)} \mathbf{Q}_{\sigma}^{j} \mathrm{d}\mathbf{x} &= \int_{V(t)} d_t \mathbf{Q}_{\sigma}^{j} + \nabla \cdot (\mathbf{Q}_{\sigma}^{j}(\mathbf{u}^{*} \mathrm{d}t + \boldsymbol{\sigma}_{t} \mathrm{d}\mathbf{B}_{t})) + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \mathbf{Q}_{\sigma}^{j}) \mathrm{d}t \mathrm{d}\mathbf{x} \nonumber \\
&= \int_{V(t)} F \mathrm{d}t \mathrm{d}\mathbf{x}. \tag{65}
\end{align}
$$

We can then write

$$
\begin{align*}
\int_{V(t)} \mathbf{Q}_{\sigma} \circ \mathrm{d}\mathbf{B}_t \, \mathrm{d}\mathbf{x} &= \int_{V(t)} \mathbf{Q}_{\sigma} \mathrm{d}\mathbf{B}_t + \frac{1}{2} d_t \left\langle \int_0^t d_t \mathbf{Q}_{\sigma}^i \cdot \int_0^t d\mathbf{B}_s^i \right\rangle \mathrm{d}\mathbf{x} \\
&= \int_{V(t)} \mathbf{Q}_{\sigma} \mathrm{d}\mathbf{B}_t - \frac{1}{2} d_t \left\langle \int_0^t \nabla \cdot (\mathbf{Q}_{\sigma}^i \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s) \cdot \int_0^t d\mathbf{B}_s^i \right\rangle \mathrm{d}\mathbf{x} \\
&= \int_{V(t)} \mathbf{Q}_{\sigma} \mathrm{d}\mathbf{B}_t - \frac{1}{2} d_t \left\langle \int_0^t \frac{\partial}{\partial x_j} (\mathbf{Q}_{\sigma}^i \boldsymbol{\sigma}_{s,jk} \mathrm{d}\mathbf{B}_s^k) \cdot \int_0^t d\mathbf{B}_s^i \right\rangle \mathrm{d}\mathbf{x} \\
&= \int_{V(t)} \mathbf{Q}_{\sigma} \mathrm{d}\mathbf{B}_t - \frac{1}{2} \frac{\partial}{\partial x_j} (\mathbf{Q}_{\sigma}^i \boldsymbol{\sigma}_{t,ji}) \mathrm{d}t \, \mathrm{d}\mathbf{x} \\
&= \int_{V(t)} \mathbf{Q}_{\sigma} \mathrm{d}\mathbf{B}_t - \frac{1}{2} \nabla \cdot (\boldsymbol{\sigma}_t \mathbf{Q}_{\sigma}) \mathrm{d}t \, \mathrm{d}\mathbf{x}. \tag{66}
\end{align*}
$$

Assembling everything and dropping the space integral, we obtain
$$
\begin{align}
& d_t \circ q + \nabla \cdot \left( q \left( \left( \boldsymbol{u} - \frac{1}{2} \nabla \cdot \boldsymbol{a} + \frac{1}{2} \boldsymbol{\sigma}_t (\nabla \cdot \boldsymbol{\sigma}_t) \right) dt + \boldsymbol{\sigma}_t \circ d\boldsymbol{B}_t \right) \right) \nonumber \\
& \qquad - Q_t \, dt - Q_\sigma \circ d\boldsymbol{B}_t \nonumber \\
& = d_t q + \nabla \cdot \left( q \left( \left( \boldsymbol{u} - \frac{1}{2} \nabla \cdot \boldsymbol{a} + \frac{1}{2} \boldsymbol{\sigma}_t (\nabla \cdot \boldsymbol{\sigma}_t) \right) dt + \boldsymbol{\sigma}_t d\boldsymbol{B}_t \right) \right) \tag{67} \\
& \qquad + \frac{1}{2} \left[ \nabla \cdot (\boldsymbol{\sigma}_t \boldsymbol{Q}_\sigma) \, dt - \nabla \cdot (\boldsymbol{a} \nabla q) \, dt - \nabla \cdot (q \boldsymbol{\sigma}_t (\nabla \cdot \boldsymbol{\sigma}_t)) \, dt \right] \nonumber \\
& \qquad - Q_t \, dt - \boldsymbol{Q}_\sigma d\boldsymbol{B}_t + \frac{1}{2} \nabla \cdot (\boldsymbol{\sigma}_t \boldsymbol{Q}_\sigma) \, dt. \nonumber
\end{align}
$$

We note that $d_t \circ q = d_t q$ since $\frac{1}{2} d_t \left( \int_0^t d_s, \int_0^t q \right) = 0$. After simplification, we obtain

$$
\begin{align}
& d_t q + \nabla \cdot \left( q \left( \left( \boldsymbol{u} - \frac{1}{2} \nabla \cdot \boldsymbol{a} \right) dt + \boldsymbol{\sigma}_t d\boldsymbol{B}_t \right) \right) + \nabla \cdot (\boldsymbol{\sigma}_t \boldsymbol{Q}_\sigma) dt \tag{68} \\
& = \frac{1}{2} \nabla \cdot (\boldsymbol{a} \nabla q) dt + Q_t dt + \boldsymbol{Q}_\sigma d\boldsymbol{B}_t, \nonumber
\end{align}
$$

which is exactly Eq. (9).

To obtain Eq. (68), we had to assume that $d\boldsymbol{Q}_\sigma$ is time-differentiable in the Lagrangian frame, which renders the demonstration slightly more restrictive concerning the shape of the forces. We do not have to perform such an assumption in Itō form, and we consider the present appendix as a sanity check of Eq. (9).

## Appendix B: Calculation Rules

### Distributivity of the Stochastic Transport Operator

The distributivity of the stochastic transport operator is detailed in this section, in the case where the stochastic transport operator is balanced by random RHS. If the evolution of two variables $f$ and $g$ are given by

$$
\begin{align}
\mathbb{D}_t f &= F_t \, dt + \boldsymbol{F}_\sigma \cdot d\boldsymbol{B}_t \nonumber \\
\mathbb{D}_t g &= G_t \, dt + \boldsymbol{G}_\sigma \cdot d\boldsymbol{B}_t, \tag{69}
\end{align}
$$

we have then

$$
\mathbb{D}_t (fg) = f \mathbb{D}_t g + g \mathbb{D}_t f + \boldsymbol{F}_\sigma \cdot \boldsymbol{G}_\sigma \, dt - (\boldsymbol{\sigma}_t (\boldsymbol{F}_\sigma) \cdot \nabla) g \, dt - (\boldsymbol{\sigma}_t (\boldsymbol{G}_\sigma) \cdot \nabla) f \, dt, \quad (70)
$$

or less formally:
$$
\begin{equation}
\begin{aligned}
\mathbb{D}_t (fg) &= f \mathbb{D}_t g + g \mathbb{D}_t f \\
&\quad + d_t \left\langle \int_0^t \mathbf{F}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t \mathbf{G}_\sigma \cdot \mathrm{d}\mathbf{B}_s \right\rangle \\
&\quad - d_t \left\langle \int_0^t \mathbf{F}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) g \right\rangle \\
&\quad - d_t \left\langle \int_0^t \mathbf{G}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) f \right\rangle.
\end{aligned}
\tag{71}
\end{equation}
$$

This relation is useful to transform the conservative form of the Navier-Stokes equations into non-conservative form.

Proof

$$
\begin{align}
\mathbb{D}_t (fg) &= \notag \\
&= f d_t g + g d_t f + d_t \langle f, g \rangle \notag \\
&\quad + f (\mathbf{u}^* \cdot \nabla) g \, dt + g (\mathbf{u}^* \cdot \nabla) f \, dt \notag \\
&\quad + \frac{1}{2} \nabla \cdot (f \mathbf{a} \nabla g + g \mathbf{a} \nabla f) \notag \\
&= f d_t g + g d_t f \notag \\
&\quad + d_t \left\langle \int_0^t - (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) f + \mathbf{F}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t - (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) g + \mathbf{G}_\sigma \cdot \mathrm{d}\mathbf{B}_s \right\rangle \notag \\
&\quad + f (\mathbf{u}^* \cdot \nabla) g \, dt + g (\mathbf{u}^* \cdot \nabla) f \, dt \notag \\
&\quad - \frac{1}{2} \left( f \nabla \cdot (\mathbf{a} \nabla g) + ((\mathbf{a} \nabla g) \cdot \nabla) f + g \nabla \cdot (\mathbf{a} \nabla f) + ((\mathbf{a} \nabla f) \cdot \nabla) g \right) \, dt. \tag{72}
\end{align}
$$

Developing only the covariation term:

$$
\begin{align}
& d_t \left\langle \int_0^t - (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) f + \mathbf{F}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t - (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) g + \mathbf{G}_\sigma \cdot \mathrm{d}\mathbf{B}_s \right\rangle = \notag \\
& \quad ((\mathbf{a}\nabla f) \cdot \nabla) g \, dt \notag \\
& \quad + d_t \left\langle \int_0^t \mathbf{F}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t \mathbf{G}_\sigma \cdot \mathrm{d}\mathbf{B}_s \right\rangle \notag \\
& \quad - d_t \left\langle \int_0^t \mathbf{F}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) g \right\rangle \notag \\
& \quad - d_t \left\langle \int_0^t \mathbf{G}_\sigma \cdot \mathrm{d}\mathbf{B}_s, \int_0^t (\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \nabla) f \right\rangle, \tag{73}
\end{align}
$$
whose right hand side is written formally

$$ (\mathbf{a}\nabla f \cdot \nabla) g \, dt + \mathbf{F}_{\sigma} \cdot \mathbf{G}_{\sigma} \, dt - (\sigma_{t}(\mathbf{F}_{\sigma}) \cdot \nabla) g \, dt - (\sigma_{t}(\mathbf{G}_{\sigma}) \cdot \nabla) f \, dt. \quad (74) $$

Substituting (74) into (72), we obtain Eq. (70).

## Work of Random Forces

For sake of clarity, we first detail calculation rules in the context of point mechanics to define the work of random forces. We consider a random force, whose impulse is a martingale $dF_{\sigma}$. We write its elementary work in a weak sense for any differentiable function $\phi(t)$ such that $\phi(0) = \phi(T) = 0$:

$$ \int_{0}^{T} \phi(t) dW_{\sigma} = \int_{0}^{T} \phi(t) \left( \frac{\partial}{\partial t} \int_{0}^{t} dF_{\sigma} \right) \cdot dX. \quad (75) $$

The random force $\frac{\partial}{\partial t} \int_{0}^{t} dF_{\sigma}$ is written here formally, since $dF_{\sigma}$ is a martingale and cannot be differentiated in time. The work in Eq. (75) can be split in two contributions: $dW_{\sigma,u}$ associated with the displacement $\mathbf{u} \, dt$, and $dW_{\sigma,\sigma}$ associated with the displacement $\sigma_{t} \, d\mathbf{B}_{t}$. We treat them separately.

$$ \begin{aligned} \int_{0}^{T} \phi(t) dW_{\sigma,u} &= \int_{0}^{T} \phi(t) \left( \frac{\partial}{\partial t} \int_{0}^{t} dF_{\sigma} \right) \cdot \mathbf{u}(x, t) \, dt. \quad (76) \\ &= - \int_{0}^{T} \left( \int_{0}^{t} dF_{\sigma} \right) \cdot (\phi(t) \mathbf{u}(x, t))' \, dt. \end{aligned} $$

The last expression is well defined and is a proper way to write this term. We can remark that $\int_{0}^{t} dF_{\sigma}$ is homogeneous to a Brownian. If we expand the following expression for a time-differentiable function $\psi(x, t)$

$$ \begin{aligned} \int_{0}^{T} d\left(\psi(x, t) \cdot \int_{0}^{t} dF_{\sigma}\right) &= \int_{0}^{T} \left(\int_{0}^{t} dF_{\sigma}\right) \cdot d\psi(x, t) + \int_{0}^{T} \psi(x, t) dF_{\sigma} \\ &= \int_{0}^{T} \left(\int_{0}^{t} dF_{\sigma}\right) \cdot \psi'(x, t) \, dt + \int_{0}^{T} \psi(x, t) dF_{\sigma}, \end{aligned} \quad (77) $$

by taking $\psi(x, t) = \phi(t)\mathbf{u}(x, t)$, we obtain
$$ 
\begin{aligned}
\int_0^T \phi(t) dW_{\sigma,u} &= \int_0^T \phi(t) u(x,t) dF_\sigma - \int_0^T d\left(\phi(t) u(x,t) \cdot \int_0^t dF_\sigma\right) \\
&= \int_0^T \phi(t) u(x,t) dF_\sigma - \phi(T) u(x,T) \cdot F_\sigma(x,T) \\
&= \int_0^T \phi(t) u(x,t) dF_\sigma.
\end{aligned} 
\quad (78) $$

We can identify

$$ dW_{\sigma,u} = u(x,t)dF_{\sigma}. \quad (79) $$

The second term $dW_{\sigma,\sigma}$ is not well defined, even in weak form. Informally, it should balance with kinetic energy of $\sigma_t dB_t$, which is not well defined (possibly infinite), and which has not been considered in the definition of total energy. Discarding this term is consistent with the derivation of the momentum in [15], where the acceleration associated with $\sigma_t dB_t$ being highly irregular is assumed to be in balance with some forces components which are equally irregular. As a consequence in our model, there is no work of the random forces associated with the Brownian motion displacement of the control surface.

## Appendix C: Displacement of a Transported Control Surface

Let us apply the SRTT to a characteristic function ($q = 1$ in $\Omega(t)$, $q = 0$ outside) transported by the flow [16], and use the divergence theorem. We obtain the volume variation associated with a control surface transported by the stochastic flow.

$$ 
\begin{aligned}
dV(t) &= d \int_{\Omega(t)} 1 \, dx = \int_{\Omega(t)} \nabla \cdot (u^* \, dt + \sigma_t dB_t) \, dx \\
&= \int_{\delta\Omega(t)} \frac{(u^* \, dt + \sigma_t dB_t) \cdot n}{dx_{d,t}} \, dS.
\end{aligned} 
\quad (80) $$

Hence, the normal displacement of the control surface is $dX_{d,t} \cdot n$, which involves the modified advection velocity. As a consequence, the modified advection velocity has to be considered for the definitions of elementary works based on surface integrals.
References

[1] H. Aluie. Scale decomposition in compressible turbulence. *Physica D: Nonlinear Phenomena*, 247(1):54–65, 2013.

[2] J. D. Anderson and J. Wendt. *Computational fluid dynamics*, volume 206. Springer, 1995.

[3] F. Auclair, L. Bordois, Y. Dossmann, T. Duhaut, A. Paci, C. Ulses, and C. Nguyen. A non-hydrostatic non-boussinesq algorithm for free-surface ocean modelling. *Ocean Modelling*, 132:12–29, 2018.

[4] W. Bauer, P. Chandramouli, B. Chapron, L. Li, and E. Mémin. Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach. *Journal of Physical Oceanography*, Feb. 2020.

[5] R. Brecht, L. Li, W. Bauer, and E. Mémin. Rotating shallow water flow under location uncertainty with a structure-preserving discretization. *Journal of Advances in Modeling Earth Systems*, 13(12):1–28, 2021.

[6] P. Chandramouli, D. Heitz, S. Laizet, and E. Mémin. Coarse large-eddy simulations in a transitional wake flow with flow models under location uncertainty. *Computers & Fluids*, 168:170–189, 2018.

[7] A. J. Chorin. A numerical method for solving incompressible viscous flow problems. *Journal of Computational Physics*, 2(1):12–26, 1967.

[8] W. K. Dewar, J. Schoonover, T. J. McDougall, and W. R. Young. Semicompressible ocean dynamics. *Journal of Physical Oceanography*, 45(1):149–156, 2015.

[9] C. Eden. Revisiting the energetics of the ocean in Boussinesq approximation. *Journal of Physical Oceanography*, 45(3):630–637, 2015.

[10] H. Giordani, R. Bourdallé-Badie, and G. Madec. An eddy-diffusivity mass-flux parameterization for modeling oceanic convection. *Journal of Advances in Modeling Earth Systems*, 12(9):e2020MS002078, 2020.

[11] D. Hoff. The zero-Mach limit of compressible flows. *Communications in mathematical physics*, 192(3):543–554, 1998.

[12] S. Kadri Harouna and E. Mémin. Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling. *Computers and Fluids*, 156:456–469, Aug. 2017.

[13] L. D. Landau and E. M. Lifshitz. *Fluid Mechanics, Course of Theoretical Physics*, volume 6. Elsevier, 2013.

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

[15] E. Mémin. Fluid flow dynamics under location uncertainty. *Geophysical & Astrophysical Fluid Dynamics*, 108(2):119–146, 2014.

[16] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part I Random transport and general models. *Geophysical and Astrophysical Fluid Dynamics*, 111(3):149–176, Apr. 2017.

[17] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading. *Geophysical and Astrophysical Fluid Dynamics*, 111(3):177–208, Apr. 2017.

[18] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions. *Geophysical and Astrophysical Fluid Dynamics*, 111(3):209–227, Apr. 2017.

[19] K. Suselj, M. J. Kurowski, and J. Teixeira. A unified eddy-diffusivity/mass-flux approach for modeling atmospheric convection. *Journal of the Atmospheric Sciences*, 76(8):2505–2537, 2019.

[20] R. Tailleux. Thermodynamics/dynamics coupling in weakly compressible turbulent stratified fluids. *International Scholarly Research Notices*, 2012, 2012.
[21] F. L. Tucciarone, E. Mémin, and L. Li. Primitive equations under location uncertainty: Analytical description and model development. *Stochastic Transport in Upper Ocean Dynamics*, page 287, 2023.

[22] G. K. Vallis. *Atmospheric and oceanic fluid dynamics*. Cambridge University Press, 2017.

[23] L. C. Woods. The thermodynamics of fluid systems. *Oxford*, 1975.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_118.png" /></figure>
# Data Driven Stochastic Primitive Equations with Dynamic Modes Decomposition

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

Check for
updates

Francesco L. Tucciarone, Etienne Mémin, and Long Li

**Abstract** As planetary flows are characterised by interaction of phenomenons in a huge range of scales, it is unaffordable today to resolve numerically the complete ocean dynamics. In this work, a stochastic version of primitive equations are implemented into the NEMO community ocean model to assess the capability of the so-called Location Uncertainty framework in representing the small scales of the ocean flows.

**Keywords** Stochastic parametrization · Ocean modelling

## 1 Introduction

Numerical resolution of planetary flows is nowadays a key tool to investigate possible future climates. The Ocean is a key actor of climate regulation and its evolution is for that reason of major interest. High resolution simulations are however extremely expensive and their usage remains limited to small domains. Large-scale simulations remain the primary tool to investigate future states of the Ocean (and of the Atmosphere as well). These simulations however do not resolve the complex interdependence of mesoscale and sub-mesoscale dynamics that characterises the global circulation, and thus great care must be put in the choice of the parametrization of all the scales that are too small to be efficiently resolved. Recent approaches incorporate noise terms to the dynamics of the flow with the goal of modelling the unresolved (and parametrised) processes, including small-scale turbulence, boundary value and scale coarsening uncertainty, as well as discretization and numerical errors. Rigorously justified methodologies have been introduced by Mémin [1] and Holm [2], providing a theoretically consistent stochastic large scale representation of the Navier-Stokes equations [3] conserving either energy

<footnote>F. L. Tucciarone (✉) · E. Mémin · L. Li
INRIA Centre de l' Université de Rennes, UMR CNRS 6625, Rennes, France
e-mail: francesco.tucciarone@inria.fr</footnote>
or circulation, respectively. Both models rely on a stochastic decomposition of the Lagrangian trajectory into a smooth-in-time component induced by the large scale velocity and a random fast-evolving uncorrelated part. From such models, a large-scale representation with a stochastic definition of the small-scale effect emerges naturally. Moreover, compared to classical large-scale deterministic modelling, the additional degree of freedom brought by the stochastic component allows us to devise new intermediate models [4, 5, 6, 7, 8]. The Location Uncertainty (LU) approach [1] has been tested within the barotropic quasi-geostrophic model, the rotating shallow water model and the surface quasi-geostrophic model, where it has proven to be more accurate in structuring the large-scale flow [4], reproducing long-term statistics [9] and providing a good trade-off between model error representation and ensemble spread [10, 11]. This work investigates the benefits of such model in the hydrostatic primitive equations, following the work of [12] in which noise based on Empirical Orthogonal Functions (EOF) was proposed. Here, a more elaborate noise defined from a Dynamics Mode Decomposition (DMD) strategy is proposed.

## 2 Location Uncertainty (LU)

The Location Uncertainty principle consists in describing the trajectory $\mathbf{X}_t$ of a fluid particle with a stochastic decomposition of the Lagrangian trajectory, represented with the following stochastic differential equation (SDE):

$$d\mathbf{X}_t = \mathbf{v}(\mathbf{X}_t, t) dt + \boldsymbol{\sigma}(\mathbf{X}_t, t) d\mathbf{B}_t, \quad (1)$$

where $\mathbf{X}: \Omega \times \mathbb{R}^+ \to \Omega$ is the fluid flow map, i.e. the trajectory followed by fluid particles starting at initial map $\mathbf{X}|_{t=0} = \mathbf{x}_0$ of the bounded domain $\Omega \subset \mathbb{R}^3$. The first component, $\mathbf{v}(\mathbf{X}_t, t)$, acts as the smooth-in-time component of the (Lagrangian) velocity of the flow, which is correlated both in space and time and associated with the integration of the equations of motion. The second component, $\boldsymbol{\sigma}(\mathbf{X}_t, t) d\mathbf{B}_t$, is a stochastic contribution (referred to as noise) that accounts for the processes that cannot be resolved at a given resolution or that have been neglected through a given numerical or physical modelling approximation. To completely define this last component, let $H$ be the Hilbert space, $H = (L^2(S), \mathbb{R}^d)$, the space of square integrable functions over $S$ with value in $\mathbb{R}^d$, with the inner product $\langle f, g \rangle_H = \int_S (f^\dagger g) dx$ and induced norm $\|f\|_H = \sqrt{\langle f, f \rangle_H}$, and let $T$ be a finite time, $T < +\infty$. In this framework $\{\mathbf{B}_t\}_{0\le t \le T}$ is defined as an $H$-valued (cylindrical) Brownian motion [13]:

$$\mathbf{B}_t = \sum_{i \in \mathbb{N}} \hat{\beta}^i \mathbf{e}_i, \quad (2)$$
where $(e_i)_{i \in \mathbb{N}}$ is a Hilbertian orthonormal basis of $H$ and $(\hat{\beta}_i)_{i \in \mathbb{N}}$ is a sequence of independent standard Brownian motions on a stochastic basis $(\Omega, \mathcal{F}, (\mathcal{F}_t)_{t \in [0, T]}, \mathbb{P})$. The noise is then properly defined as the application of an Hilbert-Schmidt symmetric integral kernel $\sigma_t f(x) = \int_S \check{\sigma}(x, y, t) f(y) dy$ to the $H$-valued cylindrical Wiener process $\mathbf{B}$ as

$$ (\sigma_t d\mathbf{B}_t)^i(x) = \int_S \check{\sigma}_{ik}(x, y, t) d\mathbf{B}_t^k(y) dy, \quad (3) $$

where the Einstein summation notation is adopted. The role of the integrable kernel $\check{\sigma}$ is to impose a fast/small scale spatial correlation. It leads to the covariance tensor $\mathbf{Q}$

$$ \begin{aligned} Q_{ij}(\mathbf{x}, \mathbf{y}, t, s) &= \mathbb{E} \left[ (\sigma_t d\mathbf{B}_t(\mathbf{x}))^i (\sigma_t d\mathbf{B}_s(\mathbf{y}))^j \right] \\ &= \delta(t-s) \int_S \check{\sigma}_{ik}(\mathbf{x}, \mathbf{z}, t) \check{\sigma}_{kj}(\mathbf{z}, \mathbf{y}, s) d\mathbf{z}, \end{aligned} $$

of the centred Gaussian process $\sigma_t d\mathbf{B}_t \sim \mathcal{N}(0, \mathbf{Q})$. The diagonal components of the covariance tensor per unit of time, defined as $\mathbf{a}(\mathbf{x}, t)\delta(t-t')dt = \mathbf{Q}(\mathbf{x}, \mathbf{x}, t, t')$, are referred to as the variance tensor. This tensor provides a measure of the strength of the noise. Notably, the variance tensor has the dimension of a viscosity in $m^2\mathbf{s}^{-1}$ and is symmetric and positive definite. The operator $\mathbf{a}$ is a compact auto-adjoint positive definite operator on $H$, that admits hence eigenfunctions $\xi_n(\cdot, t)$ with (strictly) positive eigenvalues $\lambda_n(t)$ satisfying $\sum_{n \in \mathbb{N}} \lambda_n(t) < +\infty$. As a consequence, the noise and the variance tensor $\mathbf{a}$ can be expressed (with another sequence of independent standard Brownian motions) through the spectral representation

$$ \sigma_t d\mathbf{B}_t(\mathbf{x}) = \sum_{n \in \mathbb{N}} \lambda^{1/2}(t) \xi_n(\mathbf{x}, t) d\beta_n \quad (4) $$

$$ \mathbf{a}(\mathbf{x}, t) = \sum_{n \in \mathbb{N}} \lambda(t) \xi_n(\mathbf{x}, t) \xi_n^*(\mathbf{x}, t). \quad (5) $$

This noise term is centred, however a modification can be applied in order to consider a Lagrangian displacement of the form

$$ d\mathbf{X}_t = [\mathbf{v}(\mathbf{X}_t, t) - \sigma_t Y_t(\mathbf{x})] dt + \sigma_t d\tilde{\mathbf{B}}_t(\mathbf{X}_t). \quad (6) $$

In contrast with (1), this decomposition sees the contribution of a centred Wiener process $\tilde{\mathbf{B}}_t$ that is drifted by a correlated component $\sigma_t Y_t(\mathbf{x})$. A proof of this statement can be given with Girsanov theorem, as a new probability measure $\tilde{\mathbb{P}}$ can be built in such a way that a non centred Wiener process as
$$ \tilde{\mathbf{B}}_t = \mathbf{B}_t + \int_0^t \mathbf{Y}_s \, \mathrm{d}s, \qquad (7) $$

where {$\{\mathbf{Y}_t\}_{t}$} is a random process shifting the process $\tilde{\mathbf{B}}_t$, remains centred on {$\Omega, \mathcal{F}, \tilde{\mathbf{P}}, \{\mathcal{F}_t\}_{0\le t\le T}\}$}. The new definition of the noise reads then

$$ \sigma_t \mathrm{d}\mathbf{B}_t(\mathbf{x}) = \sigma_t \mathrm{d}\tilde{\mathbf{B}}_t(\mathbf{x}) - \sigma_t Y_t(\mathbf{x}) \mathrm{d}t, \qquad (8) $$

with {$\{\tilde{\mathbf{B}}_t\}_{t}$} a $\tilde{\mathbf{P}}$-Wiener process. All the arguments provided in the following will hold for this process under $\tilde{\mathbf{P}}$, but the usage of a drifted noise is of paramount importance when the phenomenon to be modelled displays a non-zero time average and the physical processes involved cannot be regarded as completely uncorrelated, like in the case of ocean eddies and gyres. In the following, Eq. (6) will define our Lagrangian trajectory and the tilde notation will be dropped for simplicity.

### 3 Stochastic Boussinesq Equations

Within the Location Uncertainty formalism the evolution of a random tracer $q$ transported along the stochastic flow is described by the stochastic Reynolds transport theorem, introduced in [1]. The rate of change of a scalar $q$, integrated within the volume $V_t$, is described by

$$ \mathrm{d} \int_{V_t} q(\mathbf{x}, t) \, \mathrm{d}\mathbf{x} = \int_{V_t} \{D_t q + q \nabla \cdot [\mathbf{v}^* \, \mathrm{d}t + \sigma_t \mathrm{d}\mathbf{B}_t]\} (\mathbf{x}, t) \, \mathrm{d}\mathbf{x}, \qquad (9) $$

and summarised by the operator

$$ D_t q = d_t q + [\mathbf{v}^* \, \mathrm{d}t + \sigma_t \mathrm{d}\mathbf{B}_t] \cdot \nabla q - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla q) \, \mathrm{d}t. \qquad (10) $$

In this formula, the first component of the right-hand side is the increment in time at a fixed location of the process $q$, that is $d_t q = q(\mathbf{X}_t, t + \mathrm{d}t) - q(\mathbf{X}_t, t)$, playing the role of a derivative in time for a non differentiable process. Encased in the square brackets there is the stochastic advection displacement, composed of a time correlated modified advection $\mathbf{v}^*$ and a fast evolving, time uncorrelated noise $\sigma_t \mathrm{d}\mathbf{B}_t$, both advecting the scalar $q$. Under the probability measure $\tilde{\mathbf{P}}$ the velocity $\mathbf{v}^*$ is defined as

$$ \mathbf{v}^* = \mathbf{v} - \frac{1}{2} \nabla \cdot \mathbf{a} + \sigma_t^* (\nabla \cdot \sigma_t) - \sigma_t Y_t, \qquad (11) $$

where $\mathbf{v}$ is the resolved component of the velocity, $\mathbf{v}_s = \frac{1}{2} \nabla \cdot \mathbf{a}$ is the effective transport velocity resulting from the noise inhomogeneities and the last term is
the Girsanov drift due to a non centred noise. With this operator it is possible to formulate the Boussinesq equations under location uncertainty as done in [12] and reported below, split into horizontal and vertical equations using the convention $\mathbf{v} = (\mathbf{u}, w)$ and with the buoyancy defined as $b = -g \frac{\rho-\rho_0}{\rho_0}$.

**Horizontalmomentum :**

$$D_t \mathbf{u} + f \mathbf{e}_3 \times \left( \mathbf{u} \, \mathrm{d}t + \frac{1}{2} \boldsymbol{\sigma}_t \, \mathrm{d}\mathbf{B}_t^H \right) = \nabla_H \left( -p' + \frac{\nu}{3} \nabla \cdot \mathbf{v} \right) \mathrm{d}t - \nabla_H \mathrm{d}p_t^\sigma \quad (12)$$

**Verticalmomentum :**

$$D_t w = \frac{\partial}{\partial z} \left( -p' + \frac{\nu}{3} \nabla \cdot \mathbf{v} \right) \mathrm{d}t - \frac{\partial}{\partial z} \mathrm{d}p_t^\sigma + b \mathrm{d}t \quad (13)$$

**Temperatureandsalinity :**

$$D_t T = \kappa_T \Delta T \, \mathrm{d}t, \quad (14)$$

$$D_t S = \kappa_S \Delta S \, \mathrm{d}t, \quad (15)$$

**Incompressibility :**

$$\nabla \cdot [\mathbf{v} - \mathbf{v}^s] = 0, \quad \nabla \cdot \boldsymbol{\sigma}_t \, \mathrm{d}\mathbf{B}_t = 0, \quad \nabla \cdot \boldsymbol{\sigma}_t Y_t = 0 \quad (16)$$

**Equationofstate :**

$$b = b(T, S, z). \quad (17)$$

In this formulation Temperature $T$ and Salinity $S$ are introduced as active tracers and transported along the stochastic flow, thus impacting the stochastic transport of momentum through the equation of state. The term $d\boldsymbol{p}_t^\sigma$ in Eqs. (12) and (13) is a martingale term representing (under the measure $\mathbb{P}$) a zero-mean turbulent pressure related to the noise, termed *stochastic pressure*. From this starting point, Primitive equations can be achieved through a hydrostatic hypothesis on the vertical acceleration equation, resulting in the two conditions

$$-\frac{\partial p'}{\partial z} + b = 0 \quad \text{and} \quad \frac{\partial \mathrm{d}p_t^\sigma}{\partial z} = 0, \quad (18)$$

the first one being the usual hydrostatic balance, the second being the result of the uniqueness of the semi-martingale decomposition. In this work and within the scaling used in [12], the stochastic pressure is constant along depth and is supposed to be in a pure geostrophic balance [10, 14], thus not impacting the rate of change of momentum. In this hydrostatic setting the vertical component of momentum is a diagnostic variable computed through integration of the incompressibility condition (16), and the large scale pressure $p'$ is obtained through vertical integration of the buoyancy term.
# 4 Methods

The proposed method is implemented in the level-coordinate free-surface primitive equation model NEMO [15] in a double-gyre configuration consisting of a 45° degrees rotated beta plane centred at ~30°N, 3180 km long, 2120 km wide and 4 km deep, bounded by vertical walls and with a flat bottom. Seasonal winds and buoyancy changes are imposed as external forcings to induce the creation of a strong jet that separates a cold sub-polar gyre from a warm sub-tropical gyre. The complete details of this configuration are given in [16, 17], and the parameters of the simulation were chosen accordingly to the reference papers. The only change is in the values of eddy viscosity and diffusivity, enhanced of a factor five to suppress aliasing in the salinity field observed for smaller values (see Table 1 for an overview of their values). In order to assess the benefits brought by this stochastic approach, two purely deterministic simulations at different resolutions, 1/27° (R27d) and 1/3° (R3d), are compared to a stochastic simulation at 1/3° (R3LU). Each simulation consists of 10 years of data, collected every 5 days and averaged over the 5 days. The R27d simulation has been spun-up for 100 years before collecting data for the LU framework. Similarly, an additional 1/9° deterministic simulation has been spun up for 100 years in similar conditions in order to construct an initial state for the deterministic and stochastic 1/3° simulations.

## 4.1 High Resolution Data Filtering

The high resolution data used to force the low resolution stochastic model need to be filtered before being used, in order to avoid the injection of energy scales that can jeopardise the stability of the simulation. The low resolution velocity fluctuations are obtained through spatial filtering of high resolution temporal fluctuations. First, a time average is defined on the high resolution fields as

$$ \bar{\mathbf{u}}_{\text{HR}}^{t}(\mathbf{x}) = \frac{1}{T} \int_{T} \mathbf{u}_{\text{HR}}(\mathbf{x}, t) \, dt, \qquad (19) $$

<table><caption>**Table 1** Parameters of the model experiments</caption><thead><tr><th></th><th>R27d</th><th>R3d</th><th>R3LU</th></tr></thead><tbody><tr><td>Horizontal resolution</td><td>1/27° (3.9 km)</td><td>1/3° (35.3 km)</td><td>1/3° (35.3 km)</td></tr><tr><td>Horizontal grid points</td><td>540×810</td><td>60×90</td><td>60×90</td></tr><tr><td>Vertical levels</td><td>30</td><td>30</td><td>30</td></tr><tr><td>Time step</td><td>5 min</td><td>20 min</td><td>20 min</td></tr><tr><td>Eddy viscosity</td><td>-5×10<sup>-9</sup> m<sup>4</sup>s<sup>-1</sup></td><td>-5×10<sup>-12</sup> m<sup>4</sup>s<sup>-1</sup></td><td>-5×10<sup>-12</sup> m<sup>4</sup>s<sup>-1</sup></td></tr><tr><td>Eddy diffusivity</td><td>-5×10<sup>-10</sup> m<sup>4</sup>s<sup>-1</sup></td><td>1500 m<sup>2</sup>s<sup>-1</sup></td><td>1500 m<sup>2</sup>s<sup>-1</sup></td></tr></tbody></table>
so to obtain with Reynolds decomposition the high resolution fluctuations:

$$ \mathbf{u}'_{HR}(\mathbf{x}, t) = \mathbf{u}_{HR}(\mathbf{x}, t) - \bar{\mathbf{u}}_{HR}^{t}(\mathbf{x}). \qquad (20) $$

The corresponding low resolution fluctuations are obtained through a band-pass filter as $\mathbf{u}'_{LR} = (\mathcal{G}_1 - \mathcal{G}_2) \mathbf{u}'_{HR}$, where $\mathcal{G}$ represents a Gaussian filter. These filtered fields have a smaller amount of energy compared to the original snapshots, and they are re-scaled to this amount as

$$ \mathbf{u}'_{LR} = \frac{\|\mathbf{u}'_{HR}\|_2}{\left\| \left[ (\mathcal{G}_1 - \mathcal{G}_2) \mathbf{u}'_{HR} \right]_{LR}^{\downarrow} \right\|_2} \left[ (\mathcal{G}_1 - \mathcal{G}_2) \mathbf{u}'_{HR} \right]_{LR}^{\downarrow}, \qquad (21) $$

where the downward arrow represents downscaling towards low resolution. The result of this procedure sees the velocity fluctuations have the same spatial structure as before but enhanced level of energy.

## 4.2 Off-Line Noise Modelling Through DMD

Dynamical Mode Decomposition is a methodology [19] to construct a proxy linear dynamical system to describe an unknown non-linear dynamics. In this paper DMD is applied to the evolution in time of the velocity fluctuations, that is thus approximated as

$$ \mathbf{u}'(\mathbf{x}, t_{i+1}) \sim A \mathbf{u}'(\mathbf{x}, t_i). \qquad (22) $$

Such (finite dimensional) linear dynamical system is known to have a general solution:

$$ \mathbf{u}'(\mathbf{x}, t) = \sum_{m=1}^{N} b_m \exp(\mu_m t) \phi_m(\mathbf{x}), \qquad (23) $$

where $\phi_m(\mathbf{x}) \in \mathbb{C}^d$ are the eigenvectors of $A$ associated to the eigenvalues $\mu_m \in \mathbb{C}$ and $b_m \in \mathbb{C}$ are amplitudes. In particular $\mu_m = \sigma_m + i\omega_m$, the real part $\sigma_m$ is the growth rate of the mode and $\omega_m$ is the periodic frequency of the mode $m$. Since the initial data are real valued fields, the eigenvectors, eigenvalues and amplitudes will be two-by-two complex conjugate, that is $\phi_{2p} = \bar{\phi}_{2p+1}$. Following the successful proposition of [18], we split the DMD modes into correlated and uncorrelated modes in order to define the Girsanov drift through the slow component of the dynamics and the random noise through the fast component. The two sets of modes, $\mathcal{M}^u$ for the uncorrelated noise and $\mathcal{M}^c$ for the correlated part are defined as
<figure><img src="image_120.png" /><figcaption>**Fig. 1** Illustration of the selection of the DMD modes. On the left, frequencies of the modes are plotted on the unitary circle; they are coloured differently to represent their characteristic physical time scale. At this point, a threshold $\tau_c = 25d$ is chosen to differentiate the correlated from the uncorrelated modes. On the right, over violet background are plotted the correlated modes, over orange background the uncorrelated modes. The amplitude threshold for the correlated mode $C^c$ is set to zero, while for uncorrelated modes $C^u$ is set to 2. The grey dots represent the set of uncorrelated modes below this threshold, that are thus discarded</figcaption></figure>

$$ \mathcal{M}^u = \left\{ m \in [1, N] : |\mu_m| \sim 1, |\omega_m| > \frac{\pi}{\tau_c}, |b_m| \ge C^u \right\}, \quad (24) $$

$$ \mathcal{M}^c = \left\{ m \in [1, N] : |\mu_m| \sim 1, |\omega_m| \le \frac{\pi}{\tau_c}, |b_m| \ge C^c \right\}, \quad (25) $$

where $\tau_c$ is a temporal separation scale between correlated and uncorrelated (usually set to a value for which a spectral gap is observed and fixed here to twentyfive days) and $C^u$, $C^c$ are empirical cut-off of amplitudes. A visual representation of the aforementioned procedure is given in Fig. 1. As the DMD modes are not orthogonal, a scaling is applied to avoid spurious effects and to make sure that the reconstructed data corresponds to an orthogonal projection onto the subspaces spanned by the set of modes contained in $\mathcal{M}^u$ and $\mathcal{M}^c$. The procedure reads as follow:

- Construct the Gramian matrix $G$ of components $g_{m,n} = \langle \phi_m, \phi_n \rangle_H$, with $m, n \in \mathcal{M}$;

- Invert the Gramian matrix and compute the dual set of modes $\Phi^* = G^{-1}\Phi$;

- define the amplitudes as the initial state data on the dual set of modes: $\varphi_m = \langle \mathbf{u}'(\mathbf{x}, t_0), \phi_m^* \rangle_H \phi_m$.

Such procedure is applied for $\mathcal{M} = \mathcal{M}^u$ and $\mathcal{M} = \mathcal{M}^c$ separately.
## 4.3 On-Line Noise Reconstruction

Inside the NEMO core, during the simulation, the noise and Girsanov drift are defined as:

$$ \sigma_{t,\theta} d\mathbf{B}_t = \sqrt{\tau_{\theta}} \sum_{m \in \mathcal{M}^u} \exp(i\omega_m t) \varphi_m(\mathbf{x}) d\beta_m, \quad (26) $$

$$ \sigma_t \mathbf{Y}_t = \bar{\mathbf{u}}^t + \sum_{m \in \mathcal{M}^c} \exp(i\omega_m t) \varphi_m(\mathbf{x}), \quad (27) $$

with $\xi_m = \sqrt{\tau_{\theta}}\varphi_m$. In Eqs. (26) and (27) $\tau_{\theta}$ is the process decorrelation time. It is supposed to be different for each component evolving in the system, as momentum, temperature and salinity do not diffuse with the same decorrelation time. The subscript $\theta$ can thus indicate momentum $\mathbf{u}$, temperature $T$ or salinity $S$ and the corresponding noise and variance will be indicated as $\sigma_{t,\theta}d\mathbf{B}_t(\mathbf{x})$ and $a_{\theta}(\mathbf{x}, t)$. The different decorrelation times are difficult to characterize precisely (as they depends in space), but their ratio can be justified by physical reasoning. The decorrelation times are chosen in such a way that $\tau_{\mathbf{u}} = \Delta t$, $\tau_T = \frac{\kappa T}{\kappa_{\mathbf{u}}} \Delta t$ and $\tau_S = \frac{\kappa S}{\kappa_{\mathbf{u}}} \Delta t$, where $\Delta t$ is the simulation time step and $\kappa$ the molecular diffusion coefficients. Each (eigen) frequency $\omega_n$ comes in pairs and each pair of complex Brownian motion are conjugates. The real and imaginary parts of the Brownian motion are independent. As such, both the noise and Girsanov drift are real-valued fields. The variance tensor of such noise remains stationary:

$$ a_{\theta}(\mathbf{x}) = \tau_{\theta} \sum_{m \in \mathcal{M}^u} \varphi_m(\mathbf{x}) \varphi_m^{\dagger}(\mathbf{x}). \quad (28) $$

After construction with the offline data through Eqs. (26) and (27), the noise $\sigma_{t,\theta}d\mathbf{B}_t$ is constrained to live on the tangent space of the isopycnal surfaces. This procedure is operationally implemented as the application of an isopycnal projection operator $P^{\rho}$

$$ P^{\rho} = I - \frac{\nabla \rho (\nabla \rho)^{\mathrm{T}}}{|\nabla \rho|^2} \quad (29) $$

to the noise. Being the density function of temperature and salinity, $\rho = \rho(T, S, z)$, the isopycnal projection operator carries information about the current state of the simulation. The projected noise $\sigma_t d\mathbf{B}_t^\rho(\mathbf{x}) = P^\rho \sigma_t d\mathbf{B}_t(\mathbf{x})$ is thus strongly tied to the evolution of the flow density.
# 5 Results

In this work we analyse the result of a single realisation. Qualitatively speaking, a coarser resolution simulation is far from being as representative as the fine resolution simulation, as can be seen in both Figs. 2 and 3, where the leftmost panel contains results from the R27 deterministic simulation, the centre panel contains the R3 deterministic simulation and the rightmost panel contains the R3 stochastic simulation. The effect of the forcing is manifested in the R27 simulation in a jet current roughly aligned with latitude (tilted by a 45 degrees angle in the domain frame). The presence of this structure and of a secondary and smaller jet stream, is visible in the reference papers [16, 17]. This feature is absent in the R3 deterministic simulation. Figure 2 depicts the average relative vorticity $\bar{\zeta}^{10Y} = (\partial_x v - \partial_y u) / f^{10Y}$ and shows primarily the difference between a high resolution simulation and a coarse resolution simulation. In the centre panel, representing the deterministic low resolution simulation R3d, there is no sign of the characteristic jet that can be seen in the left panel as a strong contraposition of opposite sign vorticity. Viscosity on the boundary creates a sequence of alternating bands of opposite vorticity, related to jet separation problem in low resolution simulations [20]. The stochastic R3LU simulation conversely presents a much better representation of the dynamics, as a jet can now be identified clearly. The time averaged sea surface height, represented in Fig. 3 shows for the high resolution simulation the characteristic geostrophic properties of the jet stream: the northern, cold sub-polar gyre is characterised by a smaller height than the southern, warm sub-tropical gyre. This characteristic is not visible in the deterministic low resolution simulation R3d, where one can only find the effects of the boundary, while the stochastic R3LU simulation presents

<figure><img src="image_121.png" /><figcaption>**Fig. 2** 10-years averaged relative vorticity $\zeta = (\partial_x v - \partial_y u) / f$ at the surface layer of the model for deterministic high-resolution (1/27°, left), for deterministic low resolution (1/3°, middle) and for stochastic low resolution (1/3°, right)</figcaption></figure>
<figure><img src="image_122.png" /><figcaption>**Fig. 3** 10-years averaged sea surface height for deterministic high-resolution (1/27°, left), for deterministic low resolution (1/3°, middle) and for stochastic low resolution (1/3°, right)</figcaption></figure>

<figure><img src="image_123.png" /><figcaption>**Fig. 4** 10-years averaged kinetic energy for deterministic high-resolution (1/27°, left), for deterministic low resolution (1/3°, middle) and for stochastic low resolution (1/3°, right)</figcaption></figure>

a much more faded picture of this process, which is not as intense as in the R27d but definitely present. Viscosity on the left boundary provides a strong constraint to the dynamics, biasing the representation of the jet stream in both low-resolution cases. Figure 4 shows the ten years average of kinetic energy, $\overline{KE}^{10Y} = (u^2 + v^2)/2^{10Y}$. From this picture is clear that, while the stochastic simulation includes much more features when compared to its deterministic counterpart, it is still suffering the influence of the boundary, affecting the position of the jet stream separation. Figure 5 shows the vertical profiles of horizontally averaged temperature at different times for the three simulation, $\overline{T}^{x,y}(z, t) = \frac{1}{A} \int_A T(x, y, z, t) dx dy$. The deterministic high resolution profile is plotted in green, the deterministic low resolution profile is plotted in orange and the stochastic low resolution simulation is plotted in dark red. From this plot it is difficult to assess the benefits of the stochastic
<figure><img src="image_124.png" /><figcaption>Fig. 5 Vertical profile of temperature at 5 equidistant points in time. In red, the stochastic simulation, in orange the deterministic simulation, in green the reference R27 simulation</figcaption></figure>

formalism to the temperature equation, as the stochastic and deterministic low resolutions simulation are close. However, no spurious vertical mixing is observed. This corresponds to an improvement with respect to [12], and is brought by the new methodology for noise generation detailed in the previous section. In particular, constraining the noise on the isopycnal surfaces tangent planes considerably reduces the spurious mixing in POD as well. Furthermore, the projection operator has the effect of localising the effects of the noise on the jet region. From a quantitative point of view the simulations are compared using the Root Mean Square

$$ \text{RMS}(f_M) = \sqrt{\mathbb{E}[(f_M)^2]}, \qquad (30) $$

providing a measure of the energy content for the variable $f_M$, and the Root Mean Square Error

$$ \text{RMSE}(f_M) = \sqrt{\mathbb{E}[(f_{R27})_M - f_M)^2}, \qquad (31) $$
<table><caption>**Table 2** Comparison of the values obtained with the diagnostic</caption><thead><tr><th></th><th></th><th>ζ</th><th>½(u² + v²)</th><th>T</th></tr></thead><tbody><tr><td rowspan="3">RMS(f̄'f̄)</td><td>R27d</td><td>3.25e-06</td><td>18.48e-02</td><td>6.72</td></tr><tr><td>R3d</td><td>4.30e-07</td><td>1.69e-02</td><td>6.56</td></tr><tr><td>R3LU</td><td><b>1.62e-06</b></td><td><b>9.21e-02</b></td><td>6.56</td></tr><tr><td rowspan="2">RMSE(f̄'f̄)</td><td>R3d</td><td><b>3.67e-06</b></td><td>17.17e-02</td><td>0.58</td></tr><tr><td>R3LU</td><td>4.18e-06</td><td><b>16.46e-02</b></td><td><b>0.54</b></td></tr><tr><td rowspan="2">GRE</td><td>R3d</td><td>547.78</td><td>2987.56</td><td>27.46</td></tr><tr><td>R3LU</td><td><b>2.92</b></td><td><b>1.46</b></td><td><b>7.49</b></td></tr></tbody></table>

Boldface is used to highlight the best performance according to the metric.

that in turns describes the energy content of the errors. In all previous equations the expected value is considered to be taken in the volume, i.e. $E[f] = \frac{1}{V} \int_V f \, dV$. Finally, defining $\bar{f}^t = \frac{1}{T} \int_T f \, dt$ the time average and $\sigma_f^2 = \frac{1}{T} \int_T (f - \bar{f}^t)^2 \, dt$ the time variance, the Relative Gaussian Entropy [21] (GRE) at a single point

$$ \text{GRE} = \frac{1}{2} \left[ \frac{(\bar{f}_{R27}^t - \bar{f}_{M}^t)^2}{\sigma_{f,M}^2} + \frac{\sigma_{f,R27}^2}{\sigma_{f,M}^2} - 1 - \ln \left( \frac{\sigma_{f,R27}^2}{\sigma_{f,M}^2} \right) \right], \quad (32) $$

measures with a single criterion both the mean and variance reconstructions. The first term on the right-hand side of GRE represents the error in the mean weighted by the variance of the model. The remaining terms measure the error in model variability and is referred to as “dispersion”. The lower this criterion the better the reconstruction. It can be observed from (32) that this criterion is minimal if, for all points, the mean is perfectly reconstructed and if the variance of the reference equals the one of the coarse model tested. These quantitative measures have been evaluated for the mean component of three quantities, namely vorticity $\zeta$, the horizontal energy of the flow $\frac{1}{2}(u^2 + v^2)$ and temperature $T$. The values of the proposed metrics for each simulation are given in Table 2. Figure 6 provides a visual representation of the behaviour of the fluctuations around the mean states of vorticity and energy, where the fluctuations are computed in time in a Reynolds splitting fashion, $f' = f - \bar{f}^t$. The first 2 years can be considered as the time required for some adjustment of the filtered and downsampled 1/9° initial condition. It can be outlined that, while the deterministic simulation cannot sustain the initial level of variability, its stochastic counterpart shows opposite behaviour, maintaining a higher variability. Concerning the benefits of the stochastic model to the dynamical quantities $\zeta$ and $\frac{1}{2}(u^2 + v^2)$, all the presented results show that the stochastic model outperforms the deterministic simulation. The mean flow contains much more energy than the deterministic counterpart, as stated by higher values of RMS($\bar{f}^t$), the RMSE($\bar{f}^t$) of the average fields seems to be reduced, with the exception of vorticity, for which a systematic bias in the positioning of the jet stream jeopardise the computation of the RMSE. In other words the $L^2$ norm is lower for a null vorticity field than for a vorticity field exhibiting clearly a meaningful jet
<figure><img src="image_125.png" /><figcaption>**Fig. 6** Comparison between the fluctuations Root Mean Square (in space) along 10 years. In blue, the fluctuation RMS for the deterministic high resolution simulation (R27d). In green, the low resolution stochastic simulation (R3LU), in orange the low resolution deterministic simulation (R3d). The solid line shows for the three cases the values of vorticity, with corresponding scale on the left; the dotted line shows for the three cases the values of Energy, with corresponding scale on the right</figcaption></figure>

but with bias. The corresponding fluctuations show a higher RMS in time during the whole simulation, thus the stochastic simulation is energetically closer to the high resolution simulation. The distance between the stochastic simulation and the high resolution simulation, as measured by the GRE, is lower for both vorticity and energy when compared to a deterministic simulation. The effect of the stochastic parametrisation to the distribution of temperature *T* is of more difficult assessment, as the metrics show very similar behaviours. The distance measured by the GRE is lower than for the stochastic case, as much as the corresponding RMSE. However, the RMS of the mean temperature are very similar between deterministic and stochastic simulation. A comparison of the temperature fluctuations RMS in time (similar to that of Fig. 6 and not presented in this paper) confirms this similarity by showing no sensible difference between the simulations.

# 6 Conclusions

The considered stochastic model has been implemented in the ocean model NEMO. A dynamical mode decomposition based noise was considered and has been shown to be beneficial to improve the variability of the coarse resolution models and to
represent temporal statistics in a more accurate fashion. The intrinsic variability of the model has been greatly enhanced for dynamical variables as vorticity and energy, as much as the qualitative behaviour of both long time average and time-snapshots appearance. The same benefits do not seem to apply to thermodynamic quantities like temperature.

**Acknowledgments** The authors acknowledge the support of the ERC EU project 856408-STUOD.

## References

1. Mémin, E: Fluid flow dynamics under location uncertainty. Geophysical and Astrophysical Fluid Dynamics 108, 119-197 (2014).

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

3. Debusche, A., Hug, B., Mémin, E.: A consistent stochastic large-scale representation of the Navier-Stokes equations, accepted for publication in Journal of Mathematical Fluid Mechanic, 2023.

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

5. Chapron, B., Dérian, P., Mémin, E., Resseguier, V.: Large-scale flows under location uncertainty: a consistent stochastic framework. QJRMS, 144(710):251–260, 2018.

6. Cintolesi, C., Mémin, E.: Stochastic Modelling of Turbulent Flows for Numerical Simulations Fluids 5, (2020).

7. Kadri Harouna, S., Mémin, E.: Stochastic representation of the Reynolds transport theorem: Revisiting large-scale modelling. Computers and Fluids 156, 456-469 (2017).

8. Pinier, B., Mémin, E., Laizet, S., Lewandowski R.: Stochastic flow approach to model the mean velocity profile of wall-bounded flows. Phys. Rev. E, 99(6):063101, 2019.

9. Bauer, W., Chandramouli, P., Li, L., Mémin, E.: Stochastic representation of mesoscale eddy effects in coarse-resolution barotropic models Ocean Modelling, 2020, 151, pp.1-50.

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

11. Resseguier, V., Li, L., Jouan, G., Dérian, P., Mémin, E., Chapron, B.: New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics Archives of Computational Methods in Engineering, 2021, 28 (1), pp.215-261. (10.1007/s11831-020-09437-x)

12. Tucciarone F.L., Mémin, E., Li, L.: Primitive Equations Under Location Uncertainty: Analytical Description and Model Development Stochastic Transport in Upper Ocean Dynamics, Springer, 2022, pp.287-300.

13. G. D. Prato and J. Zabczyk. *Stochastic equations in infinite dimensions*. Cambridge University Press, 1992.

14. Resseguier, V., Mémin, E., Chapron, B.: Geophysical flows under location uncertainty, Part I Random transport and general models Geophysical and Astrophysical Fluid Dynamics 111, 149-176 (2017).

15. Madec, G., Bourdallé-Badie, R., Chanut, J., Clementi, E., Coward, A., Ethé, C., Iovino, D., Lea, D., Lévy, C., Lovato, T., Martin, N., Masson, S., Mocavero, S., Rousset, C., Storkey, D., Vancoppenolle, M., Müller, S., Nurser, G., Bell, M., Samson, G.: Nemo ocean engine, Oct. 2019.
16. Lévy, M., Klein, P., Tréguier, A.-M., Iovino, D., Madec, G., Masson, S., Takahashi, K.: Modifications of gyre circulation by sub-mesoscale physics. Ocean Modelling, 34(1-2):1-15, 2010.

17. Lévy, M., Resplandy, L., Klein, P., Capet, X., Iovino, D., Ethé, C.: Grid degradation of submesoscale resolving ocean models: Benefits for offline passive tracer transport. Ocean Modelling, 48:1-9, 2012.

18. Li, L., Mémin, E., Tissot, G.: Stochastic Parameterization with Dynamic Mode Decomposition Stochastic Transport in Upper Ocean Dynamics, Springer, 2022, pp.179-193.

19. Kutz J.N., Brunton S.L., Brunton, B.W., Proctor J.L.: Dynamic Mode Decomposition Society for Industrial and Applied Mathematics, 2016.

20. Chassignet, E. P. and Marshall, D. P.: Gulf Stream separation in numerical ocean models. Geophysical Monograph Series, 177, 2008.

21. Grooms, I., Majda, A., Smith, S.: Stochastic superparameterization in a quasi-geostrophic model of the antarctic circumpolar current. Ocean Modelling, 85, 10 2015.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

<figure><img src="image_126.png" /></figure>
Index

B

Bella, A., 1

C

Chapron, B., 159, 223, 247, 281
Cifani, P., 17
Collard, F., 247
Crisan, D., 53, 159

D

Dumas, F., 143

E

Ephrati, S., 17

F

Fablet, R., 247
Flandoli, F., 29

G

Garreau, P., 143
Gaultier, L., 247
Geurts, B., 193
Goodair, D., 53

H

Holm, D.D., 111, 159
Hu, R., 111

J

Jamet, Q., 143, 293

L

Lahaye, N., 1, 207, 223
Lang, O., 159
Li, L., 143, 223, 321
Lobbe, A., 159
Luesink, E., 193

M

Maingonnat, I., 207
Mémin, E., 143, 159, 223, 293, 321
Morlacchi, S., 29

O

Ouala, S., 247

P

Papini, A., 29

R

Reich, S., 261
Resseguier, V., 281
S

Street, O.D., 111

T

Tissot, G., 1, 207, 223, 293
Tucciarone, F.L., 321

V

Viviani, M., 17

Z

Zhen, Y., 281