# Data Driven Stochastic Primitive Equations with Dynamic Modes Decomposition

Francesco Tucciarone, Etienne Mémin, Long Li

► To cite this version:

Francesco Tucciarone, Etienne Mémin, Long Li. Data Driven Stochastic Primitive Equations with Dynamic Modes Decomposition. STUOD 2022 - Stochastic Transport in Upper Ocean Dynamics Annual Workshop, Sep 2022, Londres, United Kingdom. pp.321-336, 10.1007/978-3-031-40094-0_15 hal-04356242

HAL Id: hal-04356242

https://inria.hal.science/hal-04356242v1

Submitted on 15 Feb 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Data Driven Stochastic Primitive Equations with Dynamic Modes Decomposition

<figure><img src="image_2.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 = v(\mathbf{X}_t, t) dt + \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, $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, $\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 $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, $M^u$ for the uncorrelated noise and $M^c$ for the correlated part are defined as
<figure><img src="image_3.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_4.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_5.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_6.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_7.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'⟩)</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'⟩)</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 $\overline{f}^t = \frac{1}{T} \int_T f \, dt$ the time average and $\sigma_f^2 = \frac{1}{T} \int_T (f - \overline{f}^t)^2 \, dt$ the time variance, the Relative Gaussian Entropy [21] (GRE) at a single point

$$ \text{GRE} = \frac{1}{2} \left[ \frac{(\overline{f}_{R27}^t - \overline{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 - \overline{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($\overline{f}^t$), the RMSE($\overline{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_8.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_9.png" /></figure>