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

HAL
open science

# Transport noise defined from wavelet transform for model-based stochastic ocean models

Francesco Tucciarone, Long Li, Etienne Mémin, Louis Thiry

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

## To cite this version:

Francesco Tucciarone, Long Li, Etienne Mémin, Louis Thiry. Transport noise defined from wavelet transform for model-based stochastic ocean models. 2024. hal-04677272

HAL Id: hal-04677272

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

Preprint submitted on 26 Aug 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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Transport noise defined from wavelet transform
for model-based stochastic ocean models

Francesco L. Tucciarone¹, Long Li¹, Etienne Mémin¹, Louis Thiry²

Centre Inria de l'Université de Rennes – UMR CNRS 6625, Av. General Leclerc,
35042 Rennes, France

Centre Inria de Paris - CS 42112, 2 Rue Simone Iff, 75589 Paris, France

**Abstract.** The simulation of planetary flows at all the scales that have a significant impact on the climate system is unachievable with nowadays computational resources. Parametrization of the scales smaller than the simulation resolution is thus crucial to correctly resolve the ocean dynamics. In this work, a novel parametrization of the subgrid scales by means of the wavelet transform is introduced in the shallow water and primitive models within the so-called Location Uncertainty framework.

**Keywords:** Stochastic parametrization, ocean modelling

# 1 Introduction

The global climate system is strongly depending on the Ocean's state, as the interaction with the Atmosphere in the forms of mutual exchange of energy fluxes of different natures and global heat redistribution plays a crucial role in the climate regulation [1]. While observations are crucial for understanding the current state of the global ocean, numerical simulation remains the only way to forecast the system and assess future states. This is fundamental for predicting meteorological and climatological events and related hazards. Large-scale simulations of the Ocean (as well as of the Atmosphere) remain the primary tool of investigation while high resolution simulations can be obtained only for small geographical domains or short integration periods. The complex interdependence of mesoscale and sub-mesoscale dynamics, however, is lost in state-of-the-art simulations when performed at scales that are too large to capture these phenomena. Most of the modeling challenges arise from the representation of these effects in a parametrized manner [2]. A novel research trend involves incorporating perturbations and noise terms into the dynamics. The goal is to enhance the variability and parameterize sub-grid processes, turbulence, boundary value uncertainty, and account for numerical and discretization errors. Along this path, two companion methodologies have been introduced by Mémin [3] and Holm [4], providing rigorously justified methodologies to define stochastic large scales representations of the Navier-Stokes equations [5] conserving energy and circulation, respectively. These two 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 displacement noise, following ideas proposed by [6–8]. The solid theoretical background allows the definition of a large-scale representation with a stochastic component representing the subgrid contribution, introducing additional degrees of freedom to be exploited in the modelling of specific phenomena (such as large scale components [10, 11], small-scale turbulence [12, 13], boundary layer effects [14] or convection processes [15]) or to devise intermediate models [17–19, 13, 14]. The Location Uncertainty (LU) model [3] has been applied to the barotropic quasi-geostrophic model [17, 20], the baroclinic quasi-geostrophic model [27], the single-layered shallow water model [22], the surface quasi-geostrophic [21], hydrostatic primitive equations [10, 11] and recently non hydrostatic Boussinesq equations [16], proving its efficiency in structuring the large-scale flow [17], reproducing long-term statistics [20] and providing a good trade-off between model error representation and ensemble spread [21, 22]. In this work, the efficacy of a wavelet representation [23, 24] for the small scale turbulence is assessed in the context of stochastic hydrostatic primitive equations following [10, 11] and in a novel stochastic multi-layered shallow water model, based on the derivation of [22] and a modified implementation of [25].

## 2 Location uncertainty (LU)

Location Uncertainty is based on a stochastic decomposition of the Lagrangian trajectory $X_t$ of the fluid particle, so that the displacement is represented by means of the stochastic differential equation (SDE)

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

where $\mathbf{X}: S \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 $S \subset \mathbb{R}^d$ ($d=2,3$). The trajectory is thus split into a smooth-in-time (Lagrangian) velocity, $\mathbf{v}_t$, and a stochastic contribution $\boldsymbol{\sigma}_t d\mathbf{B}_t$, referred to as noise, that is non-smooth in time. The first component in equation (1) is associated to the resolved velocity in the integration of the equations of motions, while the second component accounts for processes that are either neglected or not representable at a given resolution. In order to specify the characteristic of this last (martingale) term, 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 inner product $\langle \mathbf{f}, \mathbf{g} \rangle_H = \int_S (\mathbf{f}^T \mathbf{g}) \, d\mathbf{x}$ and induced norm $\|\mathbf{f}\|_H = \langle \mathbf{f}, \mathbf{f} \rangle_H^{1/2}$ and let $T$ be a finite time, $T < +\infty$. In this context, $\{\mathbf{B}_t\}_{0 \le t \le T}$ is a cylindrical Wiener process defined on $H$ [26]:

$$\mathbf{B}_t = \sum_{i \in \mathbb{N}} \hat{\beta}^i \mathbf{e}_i, \quad (2)$$

where $(\mathbf{e}_i)_{i \in \mathbb{N}}$ is a Hilbertian orthonormal basis of the space $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 application of a Hilbert-Schmidt symmetric integral
kernel $\sigma_t f(\mathbf{x}) = \int_S \tilde{\sigma}(\mathbf{x}, \mathbf{y}, t) f(\mathbf{y}) \, d\mathbf{y}$ to the Wiener process $\mathbf{B}$ on $H$ constitutes the theoretical definition of the noise term:

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

where the Einstein convention for summation over repeated indices is adopted. The kernel $\tilde{\sigma}$ is a Hilbert-Schmidt integration kernel, assumed to be bounded in space and time. It follows that the convolution of $\tilde{\sigma}$ with $\mathbf{B}_t$ is Hilbert-Schmidt, compact, self-adjoint, positive definite and thus, by Mercer's theorem, it admits eigenfunctions and eigenvalues decreasing toward zero. This defines a centred Gaussian process

$$ \int_0^t \sigma_s \mathrm{d}\mathbf{B}_s(\mathbf{X}_s) \sim \mathcal{N}\left(0, \int_0^t \mathbf{Q}(\mathbf{X}_s, \mathbf{X}_s, s, s) \, \mathrm{d}s\right), \quad (4) $$

where the covariance tensor $\mathbf{Q}$ is defined as

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

with the integral kernel $\tilde{\sigma}$ modelled in such a way that a spatial correlation to the fast/small scale components is imposed. The strength of the noise is measured by the diagonal components of the covariance tensor per unit of time, $\mathbf{a}(\mathbf{x}, t)\delta(t-t')\mathrm{d}t = \mathbf{Q}(\mathbf{x}, \mathbf{x}, t, t')$, also referred to as the variance tensor. Notably, the variance tensor has the dimension of a viscosity in $\mathrm{m}^2\mathrm{s}^{-1}$ and is symmetric and positive definite. Furthermore, the covariance operator $\mathbf{Q}$ is a compact self-adjoint positive definite operator on $H$, that thus admits a set of orthonormal eigenfunctions $\{\xi_n(\cdot, t), n \in \mathbb{N}\}$ with (strictly) positive eigenvalues $\lambda_n(t)$ decreasing toward zero and satisfying $\sum_{n \in \mathbb{N}} \lambda_n(t) < +\infty$. Consequently, the noise term and the variance tensor can be expressed with respect to the basis provided by the eigenfunctions randomized by a series of scalar Brownian variables, $\beta_{t,n}$, as

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

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

The noise term defined above is centred, but as introduced in [27, 28] and applied in [11, 29], a modification can be applied through Girsanov transformation in order to consider a Lagrangian displacement of the form

$$ \mathrm{d}\mathbf{X}_t = [\mathbf{v}_t - \sigma_t \mathbf{Y}_t] \mathrm{d}t + \sigma_t \mathrm{d}\mathbf{B}_t, \quad (7) $$

where a correlated component $\sigma_t Y_t$ can be introduced to model phenomena displaying a non-zero time average like in the case of ocean eddies and gyres.
The transition from the Lagrangian point of view to the Eulerian point of view is provided by the stochastic Reynolds transport theorem (SRTT), introduced in [3]. It describes the rate of change of a random scalar $q$ transported by the stochastic flow (1) within a flow volume $V_t$:

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

with the operator

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

defining the stochastic transport operator. Each term of this operator has a physical interpretation. Proceeding in order, the first term of the right-hand side of (9) is the increment in time at a fixed location of the random process $q$, that is $d_t q = q(\mathbf{x}, t + d\mathbf{t}) - q(\mathbf{x}, t)$. This contribution plays the role of the partial time derivative for a process that is not time differentiable. In the square brackets it is enclosed the stochastic advection displacement. It involves a time correlated modified advection velocity,

$$\mathbf{v}^* = \mathbf{v} - \frac{1}{2} \nabla \cdot \mathbf{a} + \sigma_t^T (\nabla \cdot \sigma_t), \quad (10)$$

and a fast evolving, uncorrelated noise $\sigma_t d\mathbf{B}_t$. The advection of the process $q$ by this term leads to a multiplicative noise which is non Gaussian. This noise is referred to as transport noise in the literature. The second term in Equation (10) represents the effective transport velocity induced by statistical inhomogeneities of the noise term, and it is referred to as Itô-Stokes drift in [17]. In the following it is denoted as $\mathbf{v}^s = \frac{1}{2} \nabla \cdot \mathbf{a}$. The last term of the transport operator is a dissipation term that depicts the mixing mechanism due to the unresolved scales. In the following, the Location Uncertainty principle will be applied to a set of two-dimensional equations, the Shallow Water system, and to a set of three-dimensional equations, the Primitive Equations model. The stochastic transport operator $D_t$ has thus to be intended as built with two-dimensional differential operators in the former case, and with three dimensional differential operators in the latter.

## 3 Noise modelling with wavelets

The modelling of the noise is chosen to enhance the accuracy and the variability of a (large-scale) simulation in representing the effect of the truncated scales through random variables. Many data-driven approaches referenced previously have been proposed to that end (see for instance [10, 11, 21, 27]). Here, our goal is to propose a model-based approach for the noise definition relying only on the current state of the simulation. Opposite to data-driven technique, the noise hence depends only on the solution. It is important to outline that this does not
violate any principle of the LU derivation. Let us note however, that the noise needs to remain smooth enough in space to guarantee the existence of martin-gale solution [5]. A wavelet is a compactly supported wave-like oscillation that is localized in time [23, 24]. Wavelet processing has the characteristic of combining data processing in the time domain and in the frequency domain, with a reasonable trade-off. The forward wavelet transform decomposes the signal u from the time domain to its representation in the wavelet basis, an oscillatory waveform that reveal many signal properties and provide a sparse representation. Conversely, the inverse transform reconstructs the signal from its wavelet representation back to the time domain. The result of this operation is a set of details $\langle u, \psi_{j,k} \rangle_{L^2}$ and a large scale component $\langle u, \phi_{C,k} \rangle_{L^2}$. These fields are then randomised with a Brownian field $B_t$ defined on each point of the computational wavelet coefficients grid, so that the noise wavelet ansatz can be defined as

$$ \sigma_t d\mathbf{B}_t(\mathbf{x}) = \sum_{k=0}^{2^C-1} \langle \mathbf{u}^{(n)}, \phi_{C,k} \rangle_{L^2} d\mathbf{B}_{t,C,k} \phi_{C,k}(\mathbf{x}) \\ + \sum_{j=C}^{F} \sum_{k=0}^{2^j-1} \langle \mathbf{u}^{(n)}, \psi_{j,k} \rangle_{L^2} d\mathbf{B}_{t,j,k} \psi_{j,k}(\mathbf{x}). \quad (11) $$

In the previous equation, $F$ and $C$ are indexes that divide the details and the large scale component. The superscript (n) emphasizes that the wavelet processing is applied to the current-state n of the simulation. The first component of the noise represents the randomised large scale dynamics, and is set to zero to represent the small scale features only and perform a spatial Reynolds-like decomposition. The definition of the variance tensor can then be based on the definition of the details. Such type of noise terms can easily be shown to be well defined. They are spatially regular and their regularity is given by the choice of the wavelet basis. The wavelet transform conveys a natural multi-resolution structure to the noise as well as a natural notion of spatial scale at each level of the multi-resolution hierarchy.

## 4 Stochastic shallow water model

A sketch of the stacked shallow water system is depicted in Figure 1. In the framework of location uncertainty, the governing equations for the k-th layer ($k = 1, \dots, N$) are formulated as follows:

$$ \text{Horizontal momentum:} \qquad D_t \boldsymbol{u}_k + f \boldsymbol{u}_k^{\perp} \, dt = (-\nabla P_k + \boldsymbol{F}_k) \, dt, \qquad (12) $$

$$ \text{Mass conservation:} \qquad D_t h_k + h_k \nabla \cdot \boldsymbol{u}_k \, dt = 0, \qquad (13) $$

where $\boldsymbol{u}_k = (u_k, v_k)$ denotes the horizontal velocity with $\boldsymbol{u}_k^{\perp} = (-v_k, u_k)$, $h_k$ stands for the variable layer thickness, $f$ is the Coriolis frequency, $P_k =$
$\sum_{\ell=0}^{k-1} g'_{\ell} \eta_{\ell}$ is the Montgomery potential, $\eta_{\ell} = \eta_b + \sum_{j=\ell+1}^{N} h_j$ represents the vertical position of interface $\ell$ with $\eta_b$ indicating the position of the bottom topography, $g'_{\ell} = g(\rho_{\ell+1} - \rho_{\ell})/\rho_1$ is the reduced gravity with layer density $\rho_{\ell}$ and gravity value $g$, and $F_k = \partial_z \tau|_k \approx (\tau_{k-1} - \tau_k)/h_k$ is the vertical stress divergence. In particular, we consider only a steady surface wind stress $\tau_0$ and a linear bottom drag stress $\tau_N$. Moreover, when discussing the shallow water model the stochastic transport operator has to be understood as a two-dimensional operator. Derivation of this model can be found in [22] while a discussion of its analytical properties has been done in [31].

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Illustration of a three-layer ocean shallow water model. Each layer $k$ has a uniform density $\rho_k$ and background height $H_k$, a horizontal momentum $u_k$ and a variable thickness $h_k = \eta_{k-1} - \eta_k$, where $\eta_k$ represents the position of the interface between layers $k$ and $k+1$.</figcaption></figure>

For numerical studies, we consider a three-layer shallow water system with a steady symmetric zonal wind stress $\tau_0^x = (\tau_0/\rho_1) \cos(2\pi(y - Y/2)/Y)$, a flat bottom $\eta_b = -\sum_{k=1}^N H_k$, and a linear bottom drag $\tau_N = (\delta_{ek} f_0/2) u_N$. The common parameters for all the simulations are listed in Table 1. Time integration is performed with a third order Strong Stability Preserving Runge-Kutta (SSPRK3) method [32] for the deterministic part and a Milstein scheme (without Levy area) for the stochastic part [33, 34]. The time step is set to $0.6\Delta x/\sqrt{-g\eta_b}$ for a given grid spacing $\Delta x$. The unresolved external gravity waves are filtered using the method proposed by [35]. Advection of deterministic fluxes is performed with a fifth order Weighted Essentially Non-Oscillatory (WENO) scheme, while a second order centred scheme is applied to transport noise. The numerical implementation of this configuration follows tightly that of MQgeometry-1.0 [25].

Three simulation are performed following a spin-up run as described below: a deterministic high resolution simulation at 5km, filtered and subsampled at 10km
<table><caption>**Table 1.** Common parameters for all the models</caption><thead><tr><td>Parameters</td><td>Value</td><td>Description</td></tr></thead><tbody><tr><td>X × Y</td><td>(5120,5120) km</td><td>Domain size</td></tr><tr><td>H<sub>k</sub></td><td>(400,1100,2600) m</td><td>Mean layer thickness</td></tr><tr><td>g'<sub>k</sub></td><td>(0.025,0.0125) m s<sup>-2</sup></td><td>Reduced gravity</td></tr><tr><td>f<sub>0</sub></td><td>9.375 × 10<sup>-5</sup> s<sup>-1</sup></td><td>Mean Coriolis</td></tr><tr><td>β</td><td>1.754 × 10<sup>-11</sup> m<sup>-1</sup> s<sup>-1</sup></td><td>Coriolis gradient</td></tr><tr><td>τ<sub>0</sub></td><td>0.08 Pa</td><td>Wind stress magnitude</td></tr><tr><td>δ<sub>ek</sub></td><td>2 m</td><td>Bottom Ekman layer thickness</td></tr></tbody></table>

resolution, is taken as a reference (and thus named REF); a deterministic coarse simulation at 10km (named DET) is taken as a reference for the low resolution model. Finally, a stochastic simulation (named STO) is performed at 10km. Considering Figure 2 it is noticeable that the proposed localised basis enhances the presence of filaments and small eddies along the meandering eastward jet.

<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Comparison of instantaneous surface vorticity (top) with the zoomed version in the jet region (bottom) provided by different models at 10 km.</figcaption></figure>

This result can be further highlighted by the temporal standard deviation of the surface relative vorticity ($\omega_1 = \nabla \times u_1$), as shown in the top row of Figure 3. We observe that the STO model produces greater low-frequency variability in
the most energetic zonal jet region than the DET model at the same resolution. However, the latter allows the jet to extend further east than the former. To maintain the jet further east for the STO model, a time-correlated unresolved flow component can be added onto the uncorrelated noise through Girsanov transformation, as successfully demonstrated in our previous works [10, 11, 27, 28]. This could be performed in future work. Additionally, as illustrated in the bottom row of Figure 3, we also observe a homogenization effect of the ocean middle layer potential (PV) in the central area for both models, which corresponds well to oceanographic theory [39, 40]. Note that the PV in the middle layer is defined as $q_2 = H_2(\omega_2 + f)/h_2$, and the magnitude of the gradient of its temporal mean ($|\nabla \bar{q}_2|$) is evaluated to measure the homogeneity.

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Comparison of the (top row) temporal standard deviation of surface layer relative vorticity and (bottom row) homogenization of time-averaged potential vorticity in the middle layer, using 10-year data provided by different models (grouped by columns). The area-integrated values of $\sigma(\omega_1)/f_0$ in the most energetic zonal jet regions (highlighted by red boxes) for the DET, STO and REF models are 0.024, 0.025 and 0.032, respectively.</figcaption></figure>

We then investigate the ensemble statistical properties of the proposed stochastic model by performing 20 random realizations. Figure 4 shows the establishment of an enstrophy transfer mechanism from the large scale mean flow towards the small scale turbulent eddies. This can be seen from the top row where the progressive decrease of the ensemble average of surface vorticity is associated to an increase of its ensemble variance.
<figure><img src="image_7.png" /><figcaption>Fig. 4. Time evolution (from left to right) of ensemble mean (top) and standard deviation (bottom) of surface vorticity provided by the stochastic model at 10 km with 20 realizations.</figcaption></figure>

We next focus on the ensemble decomposition of kinetic energy (KE) and available potential energy (APE) for the random shallow water system. Recall that the KE and APE densities for the k-th layer ($k = 1, \dots, N$) and $\ell$-th interface ($\ell = 0, \dots, N-1$) are defined as follows:

$$ KE_k = \frac{1}{2} h_k |\mathbf{u}_k|^2, \quad APE_\ell = \frac{1}{2} g'_\ell \zeta_\ell^2, \quad \zeta_\ell = \sum_{j=N}^{\ell+1} (h_j - H_j), \quad (14) $$

where $\zeta_\ell$ represents the deviation of the interface. We decompose the random thickness into $h_k = \overline{h_k} + h'_k$, where $\overline{h} = \mathbb{E}[h]$ denotes the ensemble mean thickness. Consequently, $\zeta_\ell = \overline{\zeta_\ell} + \zeta'_\ell$, allowing us to define the mean potential energy (MPE) and eddy potential energy (EPE) densities as follows:

$$ MPE_\ell = \frac{1}{2} g'_\ell \overline{\zeta_\ell}^2, \quad EPE_\ell = \frac{1}{2} g'_\ell (\zeta'_\ell)^2. \quad (15) $$

Decomposing next the momentum by $\mathbf{u}_k = \hat{\mathbf{u}}_k + \mathbf{u}_k''$ with $\hat{\mathbf{u}}_k = \overline{h_k \mathbf{u}_k}/\overline{h_k}$ the thickness-weighted momentum, we define the mean kinetic energy (MKE) and eddy kinetic energy (EKE) densities by

$$ MKE_k = \frac{1}{2} \overline{h_k} |\hat{\mathbf{u}}_k|^2, \quad EKE_k = \frac{1}{2} \overline{h_k} |\mathbf{u}_k''|^2. \quad (16) $$

Figure 5 shows the behaviour in time of these energy components. Both KE and PE (MPE is not shown as it has a different order of magnitude, but follows similar profile to MKE) are first transferred from the mean to the eddy components
within the initial integration period (2 years approximately). After this balancing time, the mean and eddy components exchange energy with each other (as can be observed from their opposite phases). This phenomenon is found to be valid both locally in the jet region and globally across the entire domain.

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Ensemble energy decomposition. Left: integrated over the whole domain; Right: integrated over the jet region.</figcaption></figure>

Figure 6 shows that the ensemble generated by the proposed stochastic model covers efficiently (within a short time) the reference solution (as checked at different spatial locations), even though the ensemble forecasts and reference start from different states (which is a normal occurrence when comparing simulation of different resolution, due to the different levels of energy sustained).

<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Spread of moderate-term forecast for ensemble runs (20 members, 10 km) compared to reference run (5 km).</figcaption></figure>
# 5 Stochastic Primitive model

Within the stochastic framework of location uncertainty the Boussinesq equations can be written as

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

$$ \text{Vertical momentum:} \\ 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 (18) $$

$$ \text{Temperature and salinity:} \\ D_t T = \kappa_T \Delta T \mathrm{d}t, \quad (19) $$

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

$$ \text{Incompressibility:} \\ \nabla \cdot [\mathbf{v} - \mathbf{v}^s] = 0, \qquad \nabla \cdot \boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t = 0, \quad (21) $$

$$ \text{Equation of state:} \\ b = b(T, S, z), \quad (22) $$

with the convention $\mathbf{v} = (\mathbf{u}, w)$ and with the buoyancy defined as $b = -g \frac{\rho - \rho_0}{\rho_0}$. As opposed to the shallow water model, within the discussion of the primitive equations the stochastic transport operator has to be intended as three-dimensional. These equations where derived in [10] using asymptotic analysis starting from the stochastic Navier-Stokes of [3]. A more recent derivation starting from compressible Navier-Stokes is provided in [16]. Temperature $T$ and Salinity $S$ are considered active tracers transported by the stochastic flow, impacting the momentum equation through the (deterministic) equation of state. Consistency between the left hand side and the forcing is provided by the term $dp_t^\sigma$ in equations (17) and (18), a martingale correction corresponding to a zero-mean turbulent pressure related to the noise, termed stochastic pressure. Primitive equations are then obtained from Boussinesq equations introducing the hydrostatic hypothesis on the vertical acceleration, that provides

$$ [\boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t - \mathbf{u}^s \mathrm{d}t] \cdot \nabla w - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) \mathrm{d}t = -\frac{\partial p}{\partial z} \mathrm{d}t - \frac{\partial \mathrm{d}p_t^\sigma}{\partial z} + b \quad (23) $$

so that the pressure and stochastic pressure can be defined in relation to the vertical component of the diagnosed large scale velocity as

$$ p'(\mathbf{x}) = \int_{\eta_b}^{z} b + \mathbf{u}^s \cdot \nabla w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) \mathrm{d}\zeta, \quad (24) $$

$$ \mathrm{d}p_t^\sigma(\mathbf{x}) = \int_{\eta_b}^{z} \boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t \cdot \nabla w \mathrm{d}\zeta. \quad (25) $$

The implementation of the stochastic Primitive Equations has been done in the level-coordinate free-surface primitive equations model NEMO [36] in a wind-forced double-gyre configuration. This setting, that has already been used in
<table><caption>**Table 2.** Parameters of the model experiments.</caption><thead><tr><th></th><th>R27d</th><th>R9d</th><th>R9LU</th></tr></thead><tbody><tr><td>Horizontal resolution</td><td>1/27° (3.9 km)</td><td>~ 1/9° (11.8 km)</td><td>~ 1/9° (11.8 km)</td></tr><tr><td>Horizontal grid points</td><td>540×810</td><td>160×256</td><td>160×256</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>15 min</td><td>15 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>-9</sup> m<sup>4</sup>s<sup>-1</sup></td><td>5×10<sup>-9</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>5×10<sup>-9</sup> m<sup>4</sup>s<sup>-1</sup></td><td>5×10<sup>-9</sup> m<sup>4</sup>s<sup>-1</sup></td></tr></tbody></table>

previous works on stochastic parameterization [10, 11], consists 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 and is fully described in [37, 38]. Table 2 summarizes the physical parameters used in the simulation, in agreement to the parameters of original papers. It has to be noticed that the resolution of the R9 simulation is slightly different from that of the original paper [37] as the wavelet noise requires the domain to be a multiple of a power of 2 when an MPI z-pencil domain decomposition strategy is employed. To assess the benefits of the stochastic parametrization two purely deterministic simulations were run at two different resolutions: 1/27° (R27d), a high resolution reference, and 1/9° (R9d), the deterministic reference. These two deterministic simulations are compared to a stochastic 1/9° simulation (R9LU). The R27d simulation has been spun-up for 100 years before collecting data for the LU framework. An initial condition for R9 has been generated starting from this simulation by filtering, downsampling and running a 10 years adjustment period. Each simulation consists of 5 years of data, collected every 5 days and averaged over the 5 days.

The effect of the stochastic parametrization is assessed on the gradient of horizontal velocity. Considering Figure 4 it can be noticed that there is an increase of the small scale structures along the jet structure and in the southern gyre, where the turbulence generated by the boundary is more intense. This effect, which is a consequence of the stochastic parametrization and the associated enhancement of the small scales variability, can be assessed with the symmetric part of the velocity gradient, the strain tensor. Recall the classical decomposition of the velocity gradient as $\nabla u = \frac{1}{2}(\nabla u + (\nabla u)^T) + \frac{1}{2}(\nabla u - (\nabla u)^T)$. From the symmetric part, the normalised strain rate can be defined for the mean flow and the fluctuations respectively, as

$$ \|S\| = \frac{1}{\|u\|_{\infty}} \left( \left| \frac{\partial u}{\partial x} \right| + \left| \frac{\partial v}{\partial x} + \frac{\partial u}{\partial y} \right| + \left| \frac{\partial v}{\partial y} \right| \right), \quad (26) $$

$$ \|s\| = \frac{1}{\|u'\|_{\infty}} \left( \left| \frac{\partial u'}{\partial x} \right| + \left| \frac{\partial v'}{\partial x} + \frac{\partial u'}{\partial y} \right| + \left| \frac{\partial v'}{\partial y} \right| \right). \quad (27) $$

The integrated total strain rate $\overline{\|S\|}^t$ and turbulent strain rate $\overline{\|s\|}^t$ provide a metric to assess the effects of the parametrization along the total duration of the
<figure><img src="image_10.png" /><figcaption>Fig. 7. Snapshot of vorticity (top) and strain rate (bottom).</figcaption></figure>
<table><caption>**Table 3.** Experimental values for total strain rate and turbulent strain rate</caption><thead><tr><th></th><th>R27d</th><th>R9d</th><th>R9LU</th><th>I</th></tr></thead><tbody><tr><td>||S||<sup>t</sup></td><td>6.7295 × 10<sup>-5</sup></td><td>2.0773 × 10<sup>-5</sup></td><td>2.5314 × 10<sup>-5</sup></td><td>9.7%</td></tr><tr><td>||s||<sup>t</sup></td><td>6.7295 × 10<sup>-5</sup></td><td>1.7570 × 10<sup>-5</sup></td><td>2.2328 × 10<sup>-5</sup></td><td>9.5%</td></tr></tbody></table>

simulation. Table 3 summarizes the numerical estimation of this improvement, that is of the order of 10% when bounded by the two deterministic formulations:

$$I(f) = \frac{f_{\text{R9LU}} - f_{\text{R9d}}}{f_{\text{R27d}} - f_{\text{R9d}}}. \qquad (28)$$

The left panel of Figure 8 shows that this increase in variability is well captured by the model across a wide range of spatial scales, leading to an increase of the turbulent energy content of the flow. Additionally, the stochastic model with the proposed noise parameterization enhances the intrinsic variability of the flow at different temporal scales, as demonstrated in the right panel of Figure 8. The deterministic large-scale simulation exhibits prominent peaks at certain frequencies (around 20 and 25 days), indicative of an over-representation of certain eddies. Additionally, the inertial slope appears to be steeper, suggesting a poor representation of the eddies’ distribution within the inertial range. These aspects are clearly rectified in the stochastic simulation: the slope is significantly weaker, almost reaching the slope of the reference. Furthermore, the anomalous peaks are attenuated, indicating a more balanced distribution of eddies within the inertial range. Eddies of larger frequencies are also better represented, suggesting a more pronounced inverse cascade since no energy is injected at this scale by the noise.

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Comparison of turbulent kinetic energy spectra over (left) spatial scales and (right) temporal scales, provided by different models.</figcaption></figure>
# 6 Conclusions

A novel, wavelet based, stochastic parametrization has been implemented in two different models to test its strengths and weaknesses. The general outcome of this study is that the addition of this model-based noise term, (that depends on the current state of the simulation and not on external data), can be beneficial in facilitating the energy transfer from large scale to small scale. Both the hydrostatic primitive equations model and the shallow water model appear to support turbulent dynamics at scales smaller than those sustained by the deterministic model. This enhancement of variability is shown to be successfully exploited in ensemble-run simulations to create a larger envelope for the spread of the shallow water model. Similar ensemble experiments with the primitive equations model will be considered in the future. For future research, we plan to investigate and incorporate a time-correlated unresolved barotropic flow component by applying Girsanov transformation to the uncorrelated noise. This addition, coupled with available observational data, aims to further enhance the accuracy and reliability of the current random model.

# 7 Acknowledgements

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

# References

1. Minobe, S., Kuwano-Yoshida, A., Komori, N. et al.: Influence of the Gulf Stream on the troposphere. Nature 452, 206–209 (2008). https://doi.org/10.1038/nature06690

2. Hewitt, H.T., Roberts, M., Mathiot, P. et al.: Resolving and Parameterising the Ocean Mesoscale in Earth System Models. Curr Clim Change Rep 6, 137–152 (2020).

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

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

5. Debussche, 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.

6. Brzeźniak, Z., Capiński, M., Flandoli, F.: Stochastic partial differential equations and turbulence. Mathematical Models and Methods in Applied Sciences 1, 41-59 (1991).

7. Brzeźniak, Z., Capiński, M., Flandoli, F.: Stochastic partial differential equations and turbulence. Stochastic Analysis and Applications 10, 523-532 (1992).

8. Mikulevicius, R., Rozovskii, B.L.: Stochastic Navier-Stokes equations for turbulent flows. SIAM J. Math. Anal. 4, 1250-1310 (2004).

9. Li, L., Deremble, B., Lahaye, N., Mémin, E.: Stochastic data-driven parameterization of unresolved eddy effects in baroclinic quasi-geostrophic model. Journal of Advances in Modeling Earth Systems, 2022
10. 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.

11. Tucciarone F.L., Mémin, E., Li, L.: Data driven stochastic primitive equations with dynamic modes decomposition. In Stochastic Transport in Upper Ocean Dynamics, Springer, pp: 321-336, 2023.

12. Chandramouli, P., Mémin, E., Heitz, D.: 4D large scale variational data assimilation of a turbulent flow with a dynamics error model. Journal of Computational Physics, Volume 412,2020

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

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

15. Quentin, J., Mémin, E., Dumas, F., Li, L., Garreau, P., Toward a Stochastic Parameterization for Oceanic Deep Convection. In Stochastic Transport in Upper Ocean Dynamics II. 2023.

16. Tissot, G., Mémin, E., Jamet, Q.: Stochastic Compressible Navier-Stokes Equations Under Location Uncertainty. In Stochastic Transport in Upper Ocean Dynamics, Springer, pp: 293-319, 2023.

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

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

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

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

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

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

23. Daubechies I.: Ten Lectures on Wavelets. Society for Industrial and Applied Mathematics (SIAM), 1992.

24. Mallat S.: A Wavelet Tour of Signal Processing (Third Edition). Academic Press. 2009.

25. Thiry, L., Li, L., Roullet, G., and Mémin, E.: MQGeometry-1.0: a multi-layer quasi-geostrophic solver on non-rectangular geometries, Revised manuscript accepted for GMD, 2024, preprint available on EGUSphere, https://doi.org/10.5194/egusphere-2023-1715, 2024.

26. G. D. Prato and J. Zabczyk. *Stochastic equations in infinite dimensions*. Cambridge University Press, 1992.
27. Li, L., Mémin, E., Tissot, G.: Stochastic Parameterization with Dynamic Mode Decomposition Stochastic Transport in Upper Ocean Dynamics, Springer, 2022, pp.179-193.

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

29. Dufée, B., Mémin, and Crisan, D., Observation-Based Noise Calibration: An Efficient Dynamics for the Ensemble Kalman Filter. Stochastic Transport in Upper Ocean Dynamics 10, Springer, 43-56, 2023.

30. Borisenko, A. I., Tarapov, I. E., Silverman, R. A.: Vector and Tensor Analysis with Applications. Dover Publications (1979).

31. Lang, O., Crisan, D., Mémin, E.: Analytical Properties for a Stochastic Rotating Shallow Water Model Under Location Uncertainty. Journal of Mathematical Fluid Mechanics, (2023).

32. Gottlieb, S., Shu, C.-W., Tadmor, E., Strong stability-preserving high-order time discretization methods. SIAM Rev. 43, 89-112 (2001)

33. Fiorini, C., Boulvard, P. M., Li, L., Mémin, E. A two-step numerical scheme in time for surface quasi geostrophic equations under location uncertainty. In Stochastic Transport in Upper Ocean Dynamics, Springer, 2022.

34. Boulvard, P. M., Mémin, E.: Diagnostic of the Lévy area for geophysical flow models in view of defining high order stochastic discrete-time schemes. Foundations of Data Science. 2023.

35. Roullet, G. and Madec, G.: Salt conservation, free surface, and varying levels: A new formulation for ocean general circulation models. Journal of Geophysical Research, 2000.

36. Madec, G., Bourdallé-Badie, R.,nut, 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., Sam-
son, G.: Nemo ocean engine, Oct. 2019.

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

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

39. Rhines, P., Young, W.: Homogenization of potential vorticity in planetary gyres. Journal of Fluid Mechanics, 122:347-367, 1982.

40. Holland, W., Keffer, T., Rhines, P.: Dynamics of the oceanic general circulation: the potential vorticity field. Nature, 308:698-705, 1984.