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

HAL
open science

# Derivation of stochastic models for coastal waves

Arnaud Debussche, Étienne Mémin, Antoine Moneyron

► To cite this version:

Arnaud Debussche, Étienne Mémin, Antoine Moneyron. Derivation of stochastic models for coastal waves. Bertrand Chapron, Dan Crisan, Darryl D. Holm, Etienne Mémin, Jane-Lisa Coughlan. Stochastic Transport in Upper Ocean Dynamics III, Springer nature, pp.183-222, 2024. hal-04632538

HAL Id: hal-04632538

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

Submitted on 2 Jul 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.

Public Domain
# Derivation of stochastic models for coastal waves

Arnaud Debussche¹, Étienne Mémin², and Antoine Moneyron *,²

¹Univ Rennes, CNRS, IRMAR UMR 6625, F35000, Rennes, France.

²Univ Rennes, INRIA, IRMAR UMR 6625, F35000, Rennes, France.

July 2, 2024

## Abstract

In this paper, we consider a stochastic nonlinear formulation of classical coastal waves models under location uncertainty (LU). In the formal setting investigated here, stochastic versions of the Serre-Green-Nagdi, Boussinesq and classical shallow water wave models are obtained through an asymptotic expansion, which is similar to the one operated in the deterministic setting. However, modified advection terms emerge, together with advection noise terms. These terms are well-known features arising from the LU formalism, based on momentum conservation principle.

## Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Stochastic flow and transport</td><td>3</td></tr><tr><td>3</td><td>Non linear Shallow Water equations under location uncertainty</td><td>5</td></tr><tr><td>3.1</td><td>Deriving the stochastic Serre-Green-Naghdi system</td><td>6</td></tr><tr><td>3.1.1</td><td>Evolution equations</td><td>6</td></tr><tr><td>3.1.2</td><td>Boundary conditions</td><td>8</td></tr><tr><td>3.1.3</td><td>The wave model</td><td>10</td></tr><tr><td>3.2</td><td>Shallow water waves approximated models</td><td>13</td></tr><tr><td>3.3</td><td>Discussion on conserved quantities</td><td>14</td></tr><tr><td>4</td><td>Numerical simulations</td><td>17</td></tr><tr><td>4.1</td><td>Qualitative analysis on the effect of the noise</td><td>19</td></tr><tr><td>4.1.1</td><td>Parameters (P<sub>1</sub>)</td><td>20</td></tr><tr><td>4.1.2</td><td>Parameters (P<sub>2</sub>)</td><td>20</td></tr><tr><td>4.2</td><td>Numerical estimation of the noise-induced spreading</td><td>21</td></tr><tr><td>5</td><td>Conclusion and discussion</td><td>22</td></tr><tr><td>A</td><td>Appendix 1: Leibniz formula</td><td>23</td></tr></table>

<footnote>*Corresponding author: antoine.moneyron@inria.fr</footnote>
# 1 Introduction

The ocean surface waves constitute an essential component of ocean dynamics, as they are directly related to a strong energy exchange with the underlying current. However, the mean current and the waves follow different dynamics, operating at different times and spatial scales. Moreover, they are based on significantly different physical assumptions: surface waves rely on potential – or irrotational – flow, whereas the turbulent dynamics of the current is expressed through a non-zero vorticity. As a result, the wave-current coupling is very complex to model numerically. For this purpose, designing simplified stochastic models of surface waves would be advantageous to account for both dynamics. This can be achieved by adding as a specific noise term in the stochastic representation of ocean circulation [1, 2].

To derive such stochastic wave models, ways to adapt the classical deterministic derivation to the stochastic setting have been widely investigated. It was shown for example that deep-ocean stochastic long waves can arise from a linearised stochastic shallow water system [3] or from a stochastic Hamiltonian formulation [4]. Water waves travelling from deep water areas to shallower regions – where the water depth is much less than their wavelength – undergo significant alterations. Swift variations in height, velocity, and direction lead to substantial modifications of the free water surface profiles. Initially resembling almost perfect sine waves, these profiles evolve into an asymmetric shape. Mathematically, coastal waves are described by different dynamical models corresponding to various approximations of the irrotational Euler equations – which are averaged along the water column ultimately.

More specifically, we use the location uncertainty principle developed in [5], which is based on the addition of a noise term on the Lagrangian formulation of the displacement. In fact, this approach is a critical aspect in ocean stochastic modelling and has been widely discussed: an SDE based stochastic generalisation of the deterministic Lagrangian expression of the flow was proposed in [6, 7] for instance, from which the authors derive an Eulerian expression. Such idea also exists in turbulence modelling: for example McWilliam and Berloff proposed stochastic parametrization built upon Langevin models of turbulence in [8], devised in the wake of Kraichnan's seminal work [9].

In this paper, we aim to derive a location uncertainty interpretation of nonlinear coastal waves models, based on the modelling principle introduced in [5]. These models are naturally accompanied by some nonlinear variabilities associated with numerical or modelling uncertainties. More fundamentally, this stochastic formalism provides a way to incorporate the effect of the non-resolved vertical turbulent component. In this work, we focus on the family of models associated with the Serre-Green-Naghdi equations [10, 11], which allow capturing the non-hydrostatic phenomena related to vertical acceleration. In particular, we aim to assess the behavior of the associated numerical simulations. Our main findings are that the LU interpretations of the Saint-Venant, Boussinesq and Serre-Green-Naghdi wave models allow to break the symmetry of the wave and introduce variability compared to their deterministic counterparts. Nevertheless, they do enjoy the same conservation principles with slight conditions on the noise structure. Hence, the LU setting would allow to explore the influence of a “conservative randomness” on the classical wave model dynamics. Plus, we provide some graphical observations of this influence using numerical simulations.

The article is organised as follows: after a brief summary on the location uncertainty principle (LU), we derive a stochastic representation of the Serre-Green-Naghdi shallow water waves model following the same physical approximations used in the deterministic setting [12]. Then, we show how simpler models can be obtained from further approximations, and discuss the conservation of usual quantities such as mass, momentum and energy. The last section describes and assesses some numerical results associated with different models.
## 2 Stochastic flow and transport

Let $\mathcal{D} \subset \mathbb{R}^d$ (with $d = 2$ or 3) denote a bounded spatial domain. The LU principle relies on a stochastic formulation of the Lagrangian trajectory $(X_t)$ of the form

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

This stochastic differential equation (SDE) decomposes the flow in terms of a smooth- in-time velocity component $v$ that is both spatially and temporally correlated, and a fast unresolved flow component $\sigma dB_t$ (called noise term in the following). This latter component must be understood in the Itō sense, and is uncorrelated in time yet correlated in space. Importantly, since the Itō integral is of null expectation (as a martingale), the relation (2.1) decomposes the flow unequivocally as a large-scale displacement component and a null mean fluctuation, the mean displacement corresponding to the large-scale component expectation.

Let us now provide a precise definition of this random fluctuation component. The noise term takes values in the Hilbert space $H := (L^2(\mathcal{D}))^d$ and is defined from the Wiener process (also called cylindrical Brownian motion) [13] on a stochastic basis $(\Omega, \mathcal{F}_t, (\mathcal{F}_t)_{t \in [0,T]}, \mathbb{P})$. By definition, this stochastic basis is composed of a probability space $(\Omega, \mathbb{P})$ with filtration $(\mathcal{F}_t)_{t \in [0,T]}$ – i.e. a non decreasing family of sigma algebras evolving in time. The noise term involves independent Wiener process components $(B_t^i)_{i=1,...,d}$ defined on an orthogonal basis $(e_n)_{n \in \mathbb{N}}$ of $H$ as

$$(B_t^i)(x) = \sum_{n \in \mathbb{N}} \beta_t^n e_n^i(x), \qquad (2.2)$$

where $(\beta_t^i)_{i \in \mathbb{N}}$ is a sequence of independent one dimensional standard Brownian motions. The spatial structure of the unresolved flow component is modelled by the correlation operator, $\sigma_t$, defined as an integral operator on $H$. Let a matrix kernel $\check{\sigma} = (\check{\sigma}_{ij})_{i,j=1,...,d}$ that is bounded in space and time, then define

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

Here, we have assumed that this operator is deterministic. However, it is important to note that, if required, it could be defined as a random operator. An example of such a random correlation operator can be found in [14] where it is defined from the dynamic mode decomposition technique (DMD) [15] and a Girsanov transform. Girsanov transform enables in particular to introduce a drift term associated to a non-centred noise, through a change of probability measure. This procedure has proven particularly useful for noise calibration in data assimilation [16]. The composition of $\sigma_t[\bullet]$ and its adjoint operator $\sigma_t^*[\bullet]$ defines a compact self-adjoint positive operator, of which eigenfunctions and eigenvalues are denoted $\xi_n(\cdot, t)$ and $\lambda_n(t)$ respectively. These eigenvalues fulfil $\sum_{n \in \mathbb{N}} \lambda_n(t) < +\infty$ and decrease toward zero at infinity, and the eigenfunctions form an orthonormal basis of $H$. As such, the noise can be equivalently defined on this basis as the following spectral expression,

$$\sigma_t dB_t(x) = \sum_{n \in \mathbb{N}} \lambda_n^{1/2}(t) \xi_n(x, t) d\beta_t^n. \qquad (2.4)$$

Consequently, the noise component is a $H$-valued Gaussian process of null mean and with bounded variance – that is $E_p[\int_0^t \sigma_s dB_s] = 0, \forall t > 0$ and $E_p[\|\int_0^t \sigma_s dB_s\|_H^2] < +\infty, \forall t > 0$ – under the probability measure $\mathbb{P}$. Moreover, the auto-covariance at point $x \in \mathcal{D}$ of the unresolved flow component at each instant $t \in [0, T]$ is given by the matrix kernel of the composite operator $\sigma\sigma^*$, and denoted by $a(x, t)$, namely

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

The process $\int_0^t a(x, s) ds$ corresponds to the quadratic variation of $\int_0^t \sigma_s dB_s(x)$ [17].

The stochastic integral defining the noise could have been defined in terms of a Stratonovich integral instead
of an Itō integral [17]. For a deterministic correlation operator, **a** boils down to the local variance of the noise, due to the Itō isometry [13]. Hence, we referred to **a** as the variance tensor, although for a random correlation operator this denomination is misleading since it is a random process. Physically, the symmetric non-negative tensor **a** represents the friction coefficients of the unresolved fluid motions, and is physically homogeneous to a viscosity with unit m²/s.

Consider an extensive random tracer Θ (e.g. temperature, salinity or buoyancy) transported by the stochastic flow, fulfilling the following conservation property along the trajectories: Θ(Χ<sub>t+δt</sub>, t + δt) = Θ(Χ<sub>t</sub>, t) with δt an infinitesimal time variation. Thus, the evolution law of Θ is given by the following stochastic partial differential equation (SPDE),

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

where D<sub>t</sub> is introduced as a stochastic transport operator and d<sub>t</sub>Θ(x) := Θ(x, t + δt) - Θ(x, t) stands for the forward time-increment of Θ at a fixed point x ∈ D. This operator D<sub>t</sub> was derived in [5] using the generalized Itō formula (also called Itō-Wentzell formula in the literature) [18] and plays the role of a transport operator in a stochastic setting. Remarkably, it encompasses physically meaningful terms [17, 19]: the two first terms correspond to the classical terms of the material derivative, while the third term describes the tracer advection by the unresolved flow component. As shown in [17, 19], the resulting (non Gaussian) multiplicative noise σdB<sub>t</sub> · ∇ Θ continuously backscatters random energy to the system through the quadratic variation 1/2 (∇Θ) · (a∇Θ) of the random tracer. The last term in (2.6) represents the tracer diffusion due to the mixing of the unresolved scales. The energy loss by the diffusion term is exactly balanced with the energy brought by the noise. This pathwise balance (i.e for any realisation) leads to tracer energy conservation and highlights the parallel between the classical material derivative and the stochastic transport operator. This will be precised below.

Under specific noise and/or variance tensor definitions, the resulting diffusion [20, 5] can be connected to the additional eddy viscosity term introduced in many large-scale circulation models [21, 22]. As an additional feature of interest, this evolution law introduces an *effective* advection velocity v* in (2.6) defined for an incompressible noise term as

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

This statistical eddy-induced velocity drift captures the action of inhomogeneity of the random field on the transported tracer and the possible divergence of the unresolved flow component. It is shown in [17] that the turbophoresis term, v<sub>s</sub> = 1/2 ∇·a, can be interpreted as a generalization of the Stokes drift associated to surface wave current and that it plays a key role in the triggering of secondary circulations such as the Langmuir circulation [23, 24]. Consequently, this velocity is termed Itō-Stokes drift (ISD) in [17]. Notice that, in the modified advection (2.7), the ISD cancels out for homogeneous random fields (since then the variance tensor is constant over space).

Many useful properties of the stochastic transport operator D<sub>t</sub> have been explored by [19, 25]. In particular, if a random tracer is transported by the incompressible stochastic flow under suitable boundary conditions, then the pathwise p-th moment (p ≥ 1) of the tracer is materially and integrally invariant, namely

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

It is worth noting that the transport equation (2.6) can be written in terms of Stratonovich integral [17], as

$$ D_t \Theta = d_t \Theta + (v^* dt + \sigma \circ dB_t) \cdot \nabla \Theta = 0. \quad (2.9) $$

The Stratonovich integral has the advantage of fulfilling a “standard” chain rule, so that the notation D<sub>t</sub> for the transport operator similar to the material derivative. Itō calculus introduces second order terms – such as
the diffusion term – which becomes implicit in the Stratonovich integral. Because of this, Stratonovich noise is not anymore a martingale and is not of null expectation. Moreover, it is possible to move safely from one integral to the other under some regularity assumptions. In the following we will use the Stratonovich notation.

## 3 Non linear Shallow Water equations under location uncertainty

In this section, we derive a stochastic representation of ocean waves in the near-shore zone, focusing on the derivation of models going from the Shallow water model to the Serre-Green-Nagdi model. Our stochastic derivation is similar to the scaling procedure described in [26]. Our derivation starts from the general 3D Euler equations in the LU setting, which read

$$d_t v + \left(v - \frac{1}{2} \nabla \cdot a\right) \cdot \nabla v dt + (\sigma \circ d\mathbf{B}_t \cdot \nabla) v = -\frac{1}{\rho} \nabla (P dt + dh_t^\sigma) - g z dt, \quad (3.1)$$

$$\nabla \cdot (v - \frac{1}{2} \nabla \cdot a) = 0, \quad (3.2)$$

denoting $v = (u, w)$ the three-dimensional velocity decomposed in terms of horizontal, $u$, and vertical, $w$, components. The pressure is denoted $P$, while $\rho$ and $g$ stand for the density and the gravitational vector directed along the vertical direction, respectively.

Shallow water conditions are characterised by a water depth $h$ being much smaller than the wave length scale $L \sim 1/|k|$, where $k$ denotes the wave number. This condition is often expressed through the quantity $\beta = kh_0 \ll 1$, which informs about the predominance of dispersion. In linear theory, the amplitude $A$ of the wave is small and tends to zero when the characteristic scale tends to infinity. Another quantity is usually introduced to measure the non-linear effects: $\epsilon = A/h_0$. Shoaling processes start for $\beta \le 1$ (i.e. when wavelength and depth have the same order) and ends when the waves break at $\epsilon \ge 1$ (an illustration is provided in figure 1, which was borrowed from [27]). Thus, shoaling requires asymptotic models with short wavelength and high wave amplitude. These adimensional numbers are used to build approximated solutions of the shallow water waves system, ranging from weakly nonlinear Boussinesq models (small amplitude regime, $\beta^2 \ll 1$ and $\epsilon = O(\beta^2)$) to the nonlinear Serre-Green-Naghdi system (large amplitude regime). As it will be shown in section 3.1, within the modeling under uncertainty setting the Serre-Green-Naghdi model reads

$$\bar{D}_t^H \eta = -h(\nabla_H \cdot (\bar{u} - \frac{1}{2}\Upsilon \epsilon \bar{u}_s))dt + \Upsilon \nabla_H \cdot \overline{\sigma \circ d\mathbf{B}_t^H}, \quad (3.3a)$$

$$\bar{D}_t^H \bar{u} + \nabla_H \eta dt - \frac{1}{h} \epsilon \beta^2 \nabla_H \left( \frac{h^3}{3} (d\bar{G}) \right) = \mathcal{O}(\beta^4, \epsilon \beta^4), \quad (3.3b)$$

with

$$d\bar{G}(x, t) = \bar{D}_t^H (\nabla_H \cdot \bar{u}) - \epsilon \nabla_H \cdot \left( (\bar{u} - \frac{1}{2}\Upsilon \epsilon \bar{u}_s) dt + \Upsilon^{1/2} \overline{\sigma \circ d\mathbf{B}_t^H} \right) \nabla_H \cdot \bar{u}. \quad (3.4)$$

The notation $D_t^H f$ stands for the Stratonovich transport operator with respect to vertically averaged stochastic flow components,

$$\bar{D}_t^H f = d_t f + \epsilon \nabla_H f \cdot \bar{u}^* dt + \Upsilon^{1/2} \epsilon \nabla_H f \cdot \overline{\sigma \circ d\mathbf{B}_t^H}, \quad (3.5)$$

with the depth averaged horizontal velocity

$$\bar{u}(x, t) = \frac{1}{h} \int_{0}^{h(x, t)} \mathbf{u}(z) dz. \quad (3.6)$$

The last left-hand side term of (3.3b) can be understood as a pressure term associated to the vertical velocity component acceleration corrected by compressibility effects. This term is quite intuitively a function of the average horizontal velocity divergence. Overall, the system constitutes a stochastic version of the Serre-Green-Naghdi equations [28, 29]: compared to the original deterministic model, it includes an additional transport noise
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of shoaling processes, borrowed from [27]. As the waves approach the coasts, both the typical wavelength and the water depth decrease, the former much faster than the latter. Hence, $\beta$ decreases and $\epsilon$ increases as the waves approach the coast. Thus, shoaling processes are characterised by moderately large values of the parameters $\epsilon$ and $\beta$ – that is, away enough from the coast to avoid breaking waves, yet near enough for the bottom topography to influence the wave dynamics.</figcaption></figure>

term and the contribution of the Itō-Stokes drift induced by the inhomogeneity of the small-scale fluctuations. A graphical illustration of the bottom topography, water depth and surface deformation is provided on figure 2.

<figure><img src="image_3.png" /><figcaption>Figure 2: Illustration of the bottom topography $\eta_b$, water depth $h$ and surface deformation $\eta$.</figcaption></figure>

## 3.1 Deriving the stochastic Serre-Green-Naghdi system

In this section, we derive a stochastic representation of the Serre-Green-Naghdi model.

### 3.1.1 Evolution equations

**Scaling relations** First, we proceed to the following adimensionalisation to define the asymptotic models. Consider $x = L\tilde{x}$, $z = h_0\tilde{z}$, $v = \epsilon C_0 \tilde{v}$ with $C_0 = \sqrt{gh_0}$ the characteristic long wave velocity. Also, time is scaled as $T \sim L/C_0$. The notation $\tilde{\cdot}$ stands for adimensioned variables. We will further assume the following incompressible assumption for the noise and the ISD,

$$ \nabla \cdot \sigma \circ dB_t = 0 \quad \text{and} \quad \nabla \cdot v_s = 0. \qquad (3.7) $$

These assumptions on the noise lead to the global energy conservation – expressed as the sum of the potential and kinetic energies – of the shallow water system [30], as well as the energy conservation of transported scalar
[19, 31]. The incompressibility condition on the flow boils down to the classical divergence-free condition

$$ \nabla \cdot \boldsymbol{v} = 0. \qquad (3.8) $$

This incompressibility condition leads to the scaling $w = \epsilon\beta C_0 \tilde{w}$ for the vertical velocity. In addition, we scale the noise term and the variance tensor $\boldsymbol{a}$ by $\Upsilon^{1/2}$ and $\Upsilon$, respectively. The incompressibility condition on the noise term writes

$$ \boldsymbol{\sigma} \mathrm{d} \boldsymbol{B}_t^h = \Upsilon^{1/2} L \boldsymbol{\sigma} \mathrm{d} \tilde{\boldsymbol{B}}_t^h, \qquad (3.9) $$

with a similar scaling on the horizontal noise displacement. This provides a scaling on the vertical noise component as

$$ \boldsymbol{\sigma} \mathrm{d} \boldsymbol{B}_t^z = \Upsilon^{1/2} \epsilon \beta L \tilde{\boldsymbol{\sigma}} \mathrm{d} \tilde{\boldsymbol{B}}_t^z. \qquad (3.10) $$

The variance tensor has the dimension of a viscosity in $m^2/s$. For a deterministic correlation tensor, $\boldsymbol{\sigma}$, it is related to the variance of the flow fluctuation, $\boldsymbol{v}' = (\boldsymbol{u}', w')^T$ through

$$ \boldsymbol{a} = \mathbb{E}(\boldsymbol{v}' \boldsymbol{v}'^T) \boldsymbol{\tau}, \qquad (3.11) $$

where $\boldsymbol{\tau}$ is a decorrelation time. Now denote $\boldsymbol{a}^{HH}$, $\boldsymbol{a}^{Hz}$ and $\boldsymbol{a}^{zz}$ the horizontal (2D) matrix variance tensor, the horizontal-vertical cross vector, and the vertical variance, respectively. Thus, we get the following scaling relations,

$$ \boldsymbol{a}^{HH} = \Upsilon \epsilon^2 C_0 L \tilde{\boldsymbol{a}}^{HH} \qquad (3.12a) $$

$$ \boldsymbol{a}^{Hz} = \Upsilon \epsilon^2 \beta C_0 L \tilde{\boldsymbol{a}}^{Hz} \qquad (3.12b) $$

$$ \boldsymbol{a}^{zz} = \Upsilon \epsilon^2 \beta^2 C_0 L \tilde{\boldsymbol{a}}^{zz}. \qquad (3.12c) $$

Note that these relations could have been inferred directly from the noise scalings and relation

$$ \boldsymbol{a} \mathrm{d} t = \mathbb{E}(\boldsymbol{\sigma} \mathrm{d} \boldsymbol{B}_t (\boldsymbol{\sigma} \mathrm{d} \boldsymbol{B}_t)^T). \qquad (3.13) $$

Additionally, the ISD component scales as

$$ \frac{1}{2} \nabla \cdot \boldsymbol{a} = (\boldsymbol{u}_s, w_s) = \Upsilon \epsilon^2 C_0 (\nabla_{\tilde{H}} \cdot \tilde{\boldsymbol{a}}^{HH} + \partial_{\tilde{z}} \tilde{\boldsymbol{a}}^{Hz}, \beta \nabla_{\tilde{H}} \cdot \tilde{\boldsymbol{a}}^{Hz} + \beta \partial_{\tilde{z}} \tilde{\boldsymbol{a}}^{zz}), \qquad (3.14) $$

where $\nabla_{\tilde{H}}$ refers to the gradient with respect to the horizontal coordinates. As in the classical setting, we will assume that the large-scale component of the flow is irrotational and thus

$$ \partial_{\tilde{y}} \tilde{u} = \partial_{\tilde{x}} \tilde{v} \quad \text{and} \quad \partial_{\tilde{z}} \tilde{u} = \beta^2 \nabla_{\tilde{H}} \cdot \tilde{w}, \qquad (3.15) $$

while the adimensional continuity equation is

$$ \nabla_{\tilde{H}} \tilde{\boldsymbol{u}} = -\partial_{\tilde{z}} \tilde{w}. \qquad (3.16) $$

In the long wave assumption, the pressure is scaled by the static pressure as follows

$$ P \sim \rho g h_0 \quad \text{and thus} \quad P = \rho g h_0 \tilde{P}. \qquad (3.17) $$

The same scaling is assumed for the turbulent pressure, with the additional noise variance scaling $\Upsilon^{1/2}$,

$$ \mathrm{d}h_t^\sigma = \Upsilon^{1/2} \rho g h_0 \mathrm{d}\tilde{h}_t^\sigma. \qquad (3.18) $$

For simplicity, from now on we will drop the tilde accentuation for the adimensional variables in the following.
**Incompressibility condition** Integrating over the z axis, the incompressibility condition (3.2) gives

$$ \int_{0}^{h(x,t)} \nabla_H \cdot (u - \frac{1}{2} \Upsilon \epsilon u_s) dz = - \left( w(h(x)) - \frac{1}{2} \Upsilon \epsilon w_s(h(x)) \right). \quad (3.19) $$

Using Leibniz formula and introducing the depth averaged horizontal effective velocity $\overline{u^*}$,

$$ \overline{u^*}(x, t) = \frac{1}{h} \int_{0}^{h(x,t)} \overline{u^*}(z) dz, \quad (3.20) $$

we have

$$ \nabla_H \cdot \left( \int_{0}^{h(x,t)} \overline{u^*} dz - \overline{u^*}(h) \cdot \nabla_H h \right) = \nabla_H \cdot (\overline{u^*} h) - \overline{u^*}(h) \cdot \nabla_H h = -w^*(h). \quad (3.21) $$

In the same way, we have for the noise term

$$ \nabla_H \cdot (h \overline{\sigma d\mathbf{B}_t^H}) - \nabla_H \cdot (h \sigma d\mathbf{B}_t^H(h)) = -\sigma d\mathbf{B}_t^z(h). \quad (3.22) $$

**Adimensionned momentum equations** Gathering the previous scaled relations, we obtain the following equations for the horizontal and vertical velocity,

$$ \epsilon \int_{t}^{t} \dot{u} + \epsilon^2 \nabla_H \cdot \left( (u - \frac{1}{2} \Upsilon \epsilon u_s) u \right) dt + \epsilon^2 \partial_z \left( (w - \frac{1}{2} \Upsilon \epsilon w_s) u \right) dt + \Upsilon^{1/2} \epsilon^2 \nabla \cdot \left( \sigma \circ d\mathbf{B}_t^H u \right) + \Upsilon^{1/2} \epsilon^2 \partial_z \left( \sigma \circ d\mathbf{B}_t^z u \right) = - \nabla_H P dt - \Upsilon^{1/2} \partial_x dh_t^\sigma, \quad (3.23) $$

and

$$ \epsilon \beta^2 d_t w + \epsilon^2 \beta^2 \left( \left( u - \frac{1}{2} \Upsilon \epsilon u_s \right) \cdot \nabla_H w \right) + \epsilon^2 \beta^2 \left( \left( w - \frac{1}{2} \Upsilon \epsilon w_s \right) \partial_z w \right) dt + \epsilon^2 \beta^2 \Upsilon^{1/2} (\sigma \circ d\mathbf{B}_t^H \cdot \nabla_H w) + \epsilon^2 \beta^2 \Upsilon^{1/2} (\sigma \circ d\mathbf{B}_t^z \partial_z w) = -(\partial_z P + 1) dt - \Upsilon^{1/2} \partial_z dh_t^\sigma. \quad (3.24) $$

The vertical momentum equation can be written in a more compact form, that is

$$ \epsilon \beta^2 D_t w = -\partial_z dP - 1dt, \quad (3.25) $$

where P represents the total pressure, sum of the finite variation and martingale pressures. In this formula, 1 stands for the rescaled gravity term.

### 3.1.2 Boundary conditions

**Bottom boundary conditions** For the boundary conditions we will assume that $w(x, z, t) = \sigma \circ d\mathbf{B}_t^z = 0$ on the bottom $z = 0$. Due to this last condition we have also $a^{Hz} = a^{zz} = 0$ on the bottom. The constraint of the vertical ISD from (3.14) implies that $\partial_z a^{zz}(x, 0, t) = 0$ on the bottom. As a consequence,

$$ \nabla \cdot a(x, 0, t) = \Upsilon \epsilon^2 C_0 (\nabla_H \cdot a^{HH}(x, 0, t) + \partial_z a^{Hz}(x, 0, t), 0). $$

The divergence-free constraint of the ISD leads to following condition on the bottom,

$$ \nabla_H \cdot \nabla_H \cdot a^{HH}(x, 0, t) = 0. \quad (3.26) $$

We will also consider that the noise and the ISD have the same characteristics as the large scale velocity near the bottom: due to the large-scale component being irrotational, the following hold in the vicinity of the bottom

$$ \partial_y \sigma \circ d\mathbf{B}_t^H = \partial_x \sigma \circ d\mathbf{B}_t^H, \quad \beta^2 \nabla_H \sigma \circ d\mathbf{B}_t^z = \partial_z \sigma \circ d\mathbf{B}_t^H \text{ at } z = \delta z, \quad (3.27) $$
and

$$ \partial_x v_s = \partial_y u_s, \quad \beta^2 \nabla_H w_s = \partial_z u_s \text{ at } z = \delta z. \qquad (3.28) $$

**Free surface boundary conditions** At the free surface $z = h(x, t)$, we denote by $\eta$ the surface elevation variation and consider the scaling

$$ z = h(x, t) = h_0(1 + \epsilon\eta(x, t)), $$ 
$$ \tilde{z} = \frac{1}{h_0}h(x, t) = (1 + \epsilon\eta(x, t)). \qquad (3.29) $$

which leads to

The evolution of the *dimensional* surface elevation is given by

$$ (w - w_s)dt + \sigma \circ d\mathbf{B}_t^z = D_t \eta, \qquad (3.30) $$

where an effective vertical velocity driving the dynamics of $\eta$ appears on the left-hand side. This velocity is composed of the vertical components of the velocity, the ISD, and an additive noise variable linked with rapid vertical oscillations. As usual in long wave approximations, we assume that the whole pressure is constant at the interface, and for sake of simplicity we consider that this constant is null, that is

$$ dP_t(x, h(x, t), t) = 0. \qquad (3.31) $$

The *adimensional* form of the free surface evolution (dropping the tilde accentuation) reads,

$$ (w(h) - \frac{1}{2}\Upsilon\epsilon w_s(h))dt + \Upsilon^{1/2} \sigma \circ d\mathbf{B}_t^z(h) = d_t\eta + \epsilon \nabla_H \eta \cdot (\mathbf{u}(h) - \frac{1}{2}\Upsilon\epsilon \mathbf{u}_s(h)) dt + \Upsilon^{1/2} \epsilon \nabla_H \eta \cdot \sigma \circ d\mathbf{B}_t^H(h). \quad (3.32) $$

Considering the averaged continuity equations (3.21) we obtain

$$ d_t \eta + \epsilon \nabla_H \eta \cdot (\bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s) dt + \Upsilon^{1/2} \epsilon \nabla_H \eta \cdot \overline{\sigma \circ d\mathbf{B}_t^H} = -h(\nabla_H \cdot (\bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s)) dt + \Upsilon^{1/2} \nabla_H \cdot \overline{\sigma \circ d\mathbf{B}_t^H}, \quad (3.33) $$

which can be written more compactly as

$$ \overline{D}_t^H \eta = -h(\nabla_H \cdot (\bar{u} - \frac{1}{2}\Upsilon\epsilon\bar{u}_s))dt + \Upsilon^{1/2}\nabla_H \cdot \overline{\sigma \circ d\mathbf{B}_t^H}. \qquad (3.34) $$

Notice that we have introduced an operator $\overline{D}_t^H$, that involves only depth average velocity components. Moreover, remark that for homogeneous noise (i.e. with no statistical dependence on space location), the variance tensor is constant over space and the ISD cancels out. In such a case, the surface elevation dynamics simplifies in the more intuitive form,

$$ d_t \eta + \epsilon \nabla_H \eta \cdot \bar{u} \, dt + \Upsilon^{1/2} \epsilon \nabla_H \eta \cdot \overline{\sigma \circ d\mathbf{B}_t^H} = -h(\nabla_H \cdot \bar{u} \, dt + \Upsilon^{1/2} \nabla_H \cdot \overline{\sigma \circ d\mathbf{B}_t^H}). \qquad (3.35) $$

Additionally, in the general case, the integrated ISD can be written in terms of horizontal quantities,

$$ \begin{aligned} h \bar{u}_s &= \left( \int_0^{h(x,t)} \nabla_H \cdot \mathbf{a}^{HH} + \mathbf{a}^{Hz}(h) \right) = \nabla_H \cdot (h \bar{a}^H) - \mathbf{a}^{HH}(h) \nabla_H h + \mathbf{a}^{Hz}(h) \\ &= \nabla_H \cdot (h \bar{a}^H) - \mathbf{a}^{HH}(h) \nabla_H h + \langle \sigma \circ d\mathbf{B}_t^z(h), \sigma \circ d\mathbf{B}_t^H(h) \rangle \\ &= \nabla_H \cdot (h \bar{a}^H) - \mathbf{a}^{HH}(h) \nabla_H h - \frac{1}{2} h \nabla_H \cdot \mathbf{a}^{HH}(h). \end{aligned} \qquad (3.36) $$
### 3.1.3 The wave model

**Depth averaged horizontal momentum equation** By taking the average of the horizontal momentum equation, and using Leibniz rules in time and space for stochastic processes (see appendix A), the continuity relation and the elevation evolution together with the velocity boundary conditions, we obtain

$$ \epsilon h \, d_t \bar{u} + \epsilon^2 h \left(\bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s\right) \cdot \nabla_H \bar{u} \, dt + \Upsilon^{1/2} \epsilon^2 h \left(\overline{\sigma \circ d \bar{B}_t^H} \cdot \nabla_H \bar{u}\right) \, d\bar{u} \\ -h \epsilon^2 \nabla_H \cdot \int_0^{h(x,t)} \left( \left(\bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s\right) \bar{u} \, dt + \Upsilon^{1/2} (\sigma \circ d \bar{B}_t^H \bar{u}) - \left(\bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s\right) \bar{u} \, dt + \Upsilon^{1/2} (\overline{\sigma \circ d \bar{B}_t^H \bar{u}}) \right) \, dz = - \int_0^h \nabla_H \, dP \, dz. \quad (3.37) $$

Now we express below the depth-averaged pressure term involved in the right-hand side of this equation.

**Dynamic pressure contribution** Decomposing $dP = pdt + dh_t^\sigma$, we obtain from Leibniz formula

$$ \int_0^{h(x,t)} \nabla_H(dP)dz = \nabla_H(hd\overline{P}) - dP(h)\nabla_H h = \nabla_H(hd\overline{P}). \quad (3.38) $$

Moreover, remind that $dP$ is given by integrating the vertical momentum equation (3.25), i.e

$$ -\partial_z dP = dt + \epsilon\beta^2 D_t w. \quad (3.39) $$

Thus, integrating (3.38) vertically from depth $z$ up to the surface, one has

$$ -p(x, z, t)dt - dh_t^\sigma(x, z, t) = (z - h)dt - \epsilon\beta^2 \int_z^h D_t w(x, z', t)dz', \quad (3.40) $$

which implies

$$ -h(\bar{p}dt+\overline{dh_t^\sigma}) = -\frac{1}{2}h^2dt - \epsilon\beta^2 \int_0^h \left(\int_z^h D_t w(x, z', t)dz'\right)dz. \quad (3.41) $$

In the previous relation, the second right-hand side term can be simplified through integration by part as follows

$$ \begin{aligned} -\epsilon\beta^2 \int_0^h \left(\int_z^h D_t w dz'\right) dz &= -\epsilon\beta^2 \int_0^h z \partial_z \left(\int_0^z D_t w dz' - \int_0^h D_t w dz'\right) dz + \left[z \left(\int_0^z D_t w dz' - \int_0^h D_t w dz'\right)\right]_0^h \\ &= -\epsilon\beta^2 \int_0^h z D_t w dz. \end{aligned} \quad (3.42) $$

Then, the horizontal momentum equation reads

$$ \begin{aligned} d_t \bar{u} + \epsilon \left( \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \cdot \nabla_H \right) \bar{u} \, dt + \Upsilon^{1/2} \epsilon (\sigma \circ d \bar{B}_t^H \cdot \nabla_H) \bar{u} \, du + \nabla_H \eta \, dt + \frac{1}{h} \epsilon \beta^2 \nabla_H \int_0^h z D_t w \, dz \\ = -\frac{\epsilon}{h} \nabla_H \cdot \int_0^h \left[ \left( u - \frac{1}{2} \Upsilon \epsilon u_s \right) u - \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \bar{u} \right] \, dz - \frac{\epsilon}{h} \Upsilon^{1/2} \nabla_H \cdot \int_0^h \left[ (\sigma \circ d \bar{B}_t^H u) - (\overline{\sigma \circ d \bar{B}_t^H u}) \right] \, dz. \end{aligned} \quad (3.43) $$

So far, no approximation has been introduced in the derivation of this equation. In order to compute the different terms of this averaged equation and build simplified systems, we introduce now some asymptotic approximations. In particular, the integral terms on the right-hand side of (3.43) involve both $u$ and its vertical average value $\bar{u}$. This means the system is not closed since one cannot deduce $u$ from $\bar{u}$ solely. Therefore, one needs to investigate reliable approximations of these integral terms, as well as of the vertical acceleration term on the left-hand side of (3.43).
**Expansion of the velocities at the bottom** Since the velocity potential is harmonic and irrotational, expanding the velocity component through a Taylor expansion at $z = 0$ gives

$$ 
\begin{aligned}
\boldsymbol{u}(\boldsymbol{x}, z, t) &= \boldsymbol{u}(\boldsymbol{x}, 0, t) + \partial_z \boldsymbol{u}(\boldsymbol{x}, 0, t)z + \frac{1}{2} \partial_{zz}^2 \boldsymbol{u}(\boldsymbol{x}, 0, t)z^2 + \frac{1}{6} \partial_{zzz}^3 \boldsymbol{u}(\boldsymbol{x}, 0, t)z^3 + \mathcal{O}(\beta^4 z^4) \\
&= \boldsymbol{u}(\boldsymbol{x}, 0, t) + \beta \underbrace{\nabla_H w(\boldsymbol{x}, 0, t)}_{=0} z + \frac{1}{2} \beta \partial_z \nabla_H w(\boldsymbol{x}, 0, t)z^2 + \frac{1}{6} \beta \partial_{zz}^2 \nabla_H w(\boldsymbol{x}, 0, t)z^3 + \mathcal{O}(\beta^4 z^4) \\
&= \boldsymbol{u}(\boldsymbol{x}, 0, t) - \frac{1}{2} \beta^2 \nabla_H \nabla_H \cdot \boldsymbol{u}(\boldsymbol{x}, 0, t)z^2 - \frac{1}{6} \beta^3 \underbrace{\Delta_H \nabla_H w(\boldsymbol{x}, 0, t)}_{=0} z^3 + \mathcal{O}(\beta^4 z^4).
\end{aligned}
\quad (3.44) $$

Notice that, due to both harmonicity and irrotationality conditions, the expansion involves only even orders terms. Hence, since $z$ at most of the order of $h_0$, we obtain

$$ \boldsymbol{u}(\boldsymbol{x}, z, t) = \boldsymbol{u}(\boldsymbol{x}, 0, t) - \frac{1}{2} \beta^2 \nabla_H \nabla_H \cdot \boldsymbol{u}(\boldsymbol{x}, 0, t)z^2 + \mathcal{O}(\beta^4). \quad (3.45) $$

Averaging this equation, the bottom horizontal velocity $\boldsymbol{u}^b(\boldsymbol{x}, t) = u(\boldsymbol{x}, 0, t)$ is expressed as

$$ \boldsymbol{u}^b(\boldsymbol{x}, t) = \bar{\boldsymbol{u}}(\boldsymbol{x}, t) + \frac{1}{6} \beta^2 h^2 \nabla_H \nabla_H \cdot \boldsymbol{u}^b(\boldsymbol{x}, t) + \mathcal{O}(\beta^4). \quad (3.46) $$

Using the following expansion for $u^b$ in (3.46),

$$ \boldsymbol{u}^b = \bar{\boldsymbol{u}}(\boldsymbol{x}, t) + \epsilon \boldsymbol{x}' \cdot \nabla_H \bar{\boldsymbol{u}}, \quad (3.47) $$

we obtain

$$ \boldsymbol{u}(\boldsymbol{x}, z, t) = \bar{\boldsymbol{u}}(\boldsymbol{x}, t) + \frac{1}{6} h^2 \beta^2 \nabla_H \nabla_H \cdot \bar{\boldsymbol{u}}(\boldsymbol{x}, t) - \frac{1}{2} \beta^2 \nabla_H \nabla_H \cdot \bar{\boldsymbol{u}}(\boldsymbol{x}, t)z^2 + \mathcal{O}(\beta^4, \epsilon\beta^4). \quad (3.48) $$

Remark that this expansion brings an additional error term of order $\mathcal{O}(\epsilon\beta^4)$.

For the ISD, proceeding in the exact same way, and assuming it is irrotational at the bottom, we obtain as well

$$ 
\begin{aligned}
\boldsymbol{u}_s(\boldsymbol{x}, z, t) &= \boldsymbol{u}_s(\boldsymbol{x}, 0, t) + \partial_z \boldsymbol{u}_s(\boldsymbol{x}, 0, t)z + \frac{1}{2} \partial_{zz}^2 \boldsymbol{u}_s(\boldsymbol{x}, 0, t)z^2 + \frac{1}{6} \partial_{z^3}^3 \boldsymbol{u}_s(\boldsymbol{x}, 0, t)z^3 + \mathcal{O}(\beta^4 z^4) \\
&= \boldsymbol{u}_s(\boldsymbol{x}, 0, t) - \frac{1}{2} \beta^2 \nabla_H \nabla_H \cdot \boldsymbol{u}_s(\boldsymbol{x}, 0, t)z^2 + \mathcal{O}(\beta^4),
\end{aligned}
\quad (3.49) $$

where we used irrotationality at the bottom and incompressibility of the ISD. Averaging along depth, and replacing in the above expression the bottom ISD, we get

$$ \boldsymbol{u}_s(\boldsymbol{x}, z, t) = \bar{\boldsymbol{u}}_s(\boldsymbol{x}, t) + \mathcal{O}(\beta^2, \epsilon\beta^2). \quad (3.50) $$

Furthermore, for the vertical velocity component $w$ we have

$$ 
\begin{aligned}
w(\boldsymbol{x}, y, t) &= \partial_z w(\boldsymbol{x}, 0, t)z + \frac{1}{2} \partial_{z^2}^2 w(\boldsymbol{x}, 0, t)z^2 + \frac{1}{6} \partial_{z^3}^3 w(\boldsymbol{x}, 0, t)z^3 + \mathcal{O}(\beta^4 z^4), \\
&= -\nabla_H \cdot \boldsymbol{u}(\boldsymbol{x}, 0, t)z - \beta^2 \underbrace{\frac{1}{2} \Delta_H w(\boldsymbol{x}, 0, t)}_{=0} z^2 + \beta^2 \underbrace{\frac{1}{6} \Delta_H \nabla_H \cdot \boldsymbol{u}(\boldsymbol{x}, 0, t)z^3}_{\beta^4 z^4} + \mathcal{O}(\beta^4 z^4), \\
&= -\nabla_H \cdot \boldsymbol{u}(\boldsymbol{x}, 0, t)z + \beta^2 \underbrace{\frac{1}{6} \Delta_H \nabla_H \cdot \boldsymbol{u}(\boldsymbol{x}, 0, t)z^3}_{\beta^4 z^4} + \mathcal{O}(\beta^4 z^4).
\end{aligned}
\quad (3.51) $$

Making use of the expression of the bottom velocity, we get

$$ w(\boldsymbol{x}, y, t) = -\nabla_H \cdot \bar{\boldsymbol{u}}(\boldsymbol{x}, t)z + \mathcal{O}(\beta^2). \quad (3.52) $$
Similarly, for the vertical component of the ISD, one has

$$
\begin{align}
w_s(\mathbf{x}, y, t) &= \partial_z w_s(\mathbf{x}, 0, t)z + \frac{1}{2} \partial_{z^2}^2 w_s(\mathbf{x}, 0, t)z^2 + \frac{1}{6} \partial_{z^3}^3 w_s(\mathbf{x}, 0, t)z^3 + \mathcal{O}(\beta^4 z^4), \nonumber \\
&= -\nabla_H \cdot \mathbf{u}_s(\mathbf{x}, 0, t)z + \beta^2 \frac{1}{6} \Delta_H \nabla_H \cdot \mathbf{u}_s(\mathbf{x}, 0, t)z^3 + \mathcal{O}(\beta^4 z^4), \tag{3.53}
\end{align}
$$

then

$$
w_s(\mathbf{x}, y, t) = -\nabla_H \cdot \overline{\mathbf{u}_s}(\mathbf{x}, t)z + \mathcal{O}(\beta^2). \quad (3.54)
$$

For the noise term, following the same procedure, the horizontal component is expressed as

$$
\sigma \circ d \mathbf{B}_t^H(\mathbf{x}, z, t) = \overline{\sigma \circ d \mathbf{B}_t^H}(\mathbf{x}, z) + \frac{1}{6} \beta^2 h^2 \nabla_H \cdot \overline{\nabla_H \cdot \sigma \circ d \mathbf{B}_t^H}(\mathbf{x}, z) - \frac{1}{2} \beta^2 \nabla_H \cdot \overline{\nabla_H \cdot \sigma \circ d \mathbf{B}_t^H}(\mathbf{x}, z) z^2 + \mathcal{O}(\beta^4, \epsilon\beta^4). \quad (3.55)
$$

and the vertical one as

$$
\sigma \circ d \mathbf{B}_t^z(\mathbf{x}, y, t) = - \nabla_H \cdot \overline{\sigma \circ d \mathbf{B}_t^H}(\mathbf{x}, t) z + \mathcal{O}(\beta^2). \quad (3.56)
$$

From these formulae, we can approximate the variance tensor components $a^H$, $a^{Hz}$ and $a^{zz}$:

$$
\begin{align}
a^H dt &= \mathbb{E}(\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H \otimes \boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H) = \mathbb{E}(\overline{\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H} \otimes \overline{\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H}) + \mathcal{O}(\beta^2, \epsilon\beta^2) \nonumber \\
&= \bar{\mathbf{a}}^H + \mathcal{O}(\beta^2, \epsilon\beta^2), \tag{3.57a}
\end{align}
$$

$$
\begin{align}
a^{Hz} dt &= \mathbb{E}(\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H \boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^z) = -\beta \mathbb{E}(\overline{\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H} \nabla_H \cdot \overline{\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H}) z + \mathcal{O}(\beta^2, \epsilon\beta^2) \nonumber \\
&= -\frac{1}{2} z \nabla_H \cdot \bar{\mathbf{a}}^H + \mathcal{O}(\beta^2, \epsilon\beta^2), \tag{3.57b}
\end{align}
$$

$$
\begin{align}
a^{zz} dt &= \mathbb{E}(\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^z \boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^z) = \mathbb{E}(\nabla_H \cdot \overline{\boldsymbol{\sigma} \mathrm{d}\mathbf{B}_t^H} \nabla_H \cdot \overline{\boldsymbol{\sigma} \circ \mathrm{d}\mathbf{B}_t^H}) z^2 + \mathcal{O}(\beta^4, \epsilon\beta^4) \nonumber \\
&= z^2 \bar{\mathbf{a}}^{\mathrm{div}} + \mathcal{O}(\beta^4), \tag{3.57c}
\end{align}
$$

With the above approximations, we observe that all the quadratic integrals scale as

$$
\int_0^h \left[ \left( u - \frac{1}{2} \Upsilon \epsilon u_s \right) u - \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \bar{u} \right] dz \sim \mathcal{O}(\beta^4, \epsilon\beta^4), \quad (3.58)
$$

$$
\int_0^h \left[ (\sigma \circ d \mathbf{B}_t^H u) - (\overline{\sigma \circ d \mathbf{B}_t^H} \overline{u}) \right] dz \sim \mathcal{O}(\beta^4, \epsilon\beta^4) \quad (3.59)
$$

Besides, the ISD (3.36) reads now

$$
\begin{align}
\mathbf{u}_s &= \frac{1}{h} (\nabla_H \cdot (h \bar{\mathbf{a}}^H) - \mathbf{a}^{HH}(h) \nabla_H h - \frac{1}{2} h \nabla_H \cdot \mathbf{a}^{HH}(h)) \nonumber \\
&= \frac{1}{2} (\nabla_H \cdot \bar{\mathbf{a}}^H) + \mathcal{O}(\beta^2, \epsilon\beta^2), \tag{3.60}
\end{align}
$$

$$
\begin{align}
w_s &= \nabla_H \cdot \mathbf{a}^{Hz} + \partial_z a^{zz} \nonumber \\
    &= -z \left( \frac{1}{2} \nabla_H \cdot \nabla_H \cdot \bar{\mathbf{a}}^H - 2 \bar{\mathbf{a}}^{\text{div}} \right) + \mathcal{O}(\beta^2, \epsilon\beta^2). \tag{3.61}
\end{align}
$$
Using the approximation of the vertical velocity (3.52), the vertical acceleration relation (3.25) reads

$$D_t w = -z \left( d_t \nabla_H \cdot \bar{u} + \epsilon \left( \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \cdot \nabla_H \right) \nabla_H \cdot \bar{u} dt - \epsilon \left( \nabla_H \cdot \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \nabla_H \cdot \bar{u} \right) dt + \epsilon \Upsilon^{1/2} \overline{\left( \sigma \circ d \bar{B}_t^H \cdot \nabla_H \right)} \nabla_H \cdot \bar{u} - \epsilon \Upsilon^{1/2} \overline{\left( \nabla_H \cdot \sigma \circ d \bar{B}_t^H \nabla_H \cdot \bar{u} \right)} \right) + \mathcal{O}(\beta^2, \epsilon\beta^2). \quad (3.62)$$

Eventually, we find the evolution equations for the surface elevation and the depth averaged velocity,

$$\overline{D_t^H} \eta = -h \left( \nabla_H \cdot \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) dt + \Upsilon \nabla_H \cdot \overline{\sigma \circ d \bar{B}_t^H} \right), \quad (3.63a)$$

$$d_t \bar{u} + \epsilon \left( \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \cdot \nabla_H \right) \bar{u} dt + \Upsilon^{1/2} \epsilon \overline{\left( \sigma \circ d \bar{B}_t^H \cdot \nabla_H \right)} \bar{u} + \nabla_H \eta dt - \frac{1}{h} \epsilon \beta^2 \nabla_H \left( \frac{h^3}{3} (d\bar{G}) \right) = \mathcal{O}(\beta^4, \epsilon\beta^4) \quad (3.63b)$$

with

$$d\overline{G}(x, t) = \left( d_t \nabla_H \cdot \bar{u} + \epsilon \left( \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \cdot \nabla_H \right) \nabla_H \cdot \bar{u} dt - \epsilon \left( \nabla_H \cdot \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \nabla_H \cdot \bar{u} \right) dt + \epsilon \Upsilon^{1/2} \overline{\left( \sigma \circ d \bar{B}_t^H \cdot \nabla_H \right)} \nabla_H \cdot \bar{u} dt - \epsilon \Upsilon^{1/2} \overline{\left( \nabla_H \cdot \sigma \circ d \bar{B}_t^H \nabla_H \cdot \bar{u} \right)} \right), \quad (3.64)$$

which can be more compactly written as

$$d\overline{G}(x, t) = \overline{D_t^H} (\nabla_H \cdot \bar{u}) - \epsilon \nabla_H \cdot \left( \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) dt + \Upsilon^{1/2} \overline{\sigma \circ d \bar{B}_t^H} \right) \nabla_H \cdot \bar{u}. \quad (3.65)$$

This last expression can readily be understood as the acceleration of the average horizontal velocity divergence corrected by compressibility effects. This system constitutes a stochastic version of the Serre-Green-Naghdi equations [28, 29]. Compared to the original deterministic model, additional noise terms are involved. Those terms correspond to small scale advection and are accompanied with a corresponding ISD correction term in the large scale advection.

Let us now exhibit some simplified models that arise by neglecting higher order terms.

## 3.2 Shallow water waves approximated models

In this section we present two stochastic representations of classical approximations of the Serre-Green-Naghdi equations, namely the Shallow Water and the Boussinesq wave models. In addition, we briefly mention a stochastic version of the Kordeveg-De Vries equation. Through the section, we will make use of the Stokes number (also called Ursel number), defined as $S = \epsilon/\beta^2$.

**Shallow water (or Saint-Venant) long waves approximation** For long waves regimes such as tidal waves, we have $\beta \ll \epsilon = A/h \ll 1$, which corresponds to a Stokes number $S \gg 1$. Neglecting in system (3.63) the terms of order higher than $\epsilon$ gives

$$\overline{D_t^H} \eta = -h \left( \nabla_H \cdot \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) dt + \Upsilon^{1/2} \nabla_H \cdot \overline{\sigma \circ d \bar{B}_t^H} \right), \quad (3.66a)$$

$$d_t \bar{u} + \epsilon \left( \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \cdot \nabla_H \right) \bar{u} dt + \Upsilon^{1/2} \epsilon \overline{\left( \sigma \circ d \bar{B}_t^H \cdot \nabla_H \right)} \bar{u} + \nabla_H \eta dt = 0. \quad (3.66b)$$

This system corresponds to a stochastic version of the 2D Shallow water model [30]. We note that the noise term is kept assuming $\Upsilon$ is of order 1 or higher. The diffusion term is in balance with the energy brought by the noise and must be kept to ensure energy conservation. For lower noise amplitude, the system boils down to the classical deterministic system. Notice that this system results from neglecting the vertical acceleration, which corresponds to the usual hydrostatic assumption. However, differently from the deterministic (linear) shallow
water system, this corresponding linear stochastic system admits dispersive waves as solutions due to the noise term [3].

**Boussinesq approximation** Assuming that $S \approx 1$ and $\beta \ll 1$, we retain only the terms of order $\epsilon$ and $\beta^2$ in the system (3.63). Thus, we obtain a stochastic interpretation of Boussinesq wave equation. This system reads

$$D_t^H \eta = -h(\nabla_H \cdot (\bar{u} - \frac{1}{2}\Upsilon\epsilon\bar{u}_s))dt + \Upsilon^{1/2}\nabla_H \cdot \overline{\sigma \circ d\mathbf{B}_t^H}, \quad (3.67a)$$

$$d_t \bar{u} + \epsilon \left( \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \cdot \nabla_H \right) \bar{u} dt + \Upsilon^{1/2} \epsilon \left( \overline{\sigma \circ d\mathbf{B}_t^H} \cdot \nabla_H \right) \bar{u} + \nabla_H \eta dt - \epsilon \beta^2 \left( \frac{h^2}{3} \nabla_H \cdot \nabla_H \right) d_t \bar{u} = \mathcal{O}(\beta^4, \epsilon\beta^4) \quad (3.67b)$$

Remark that an additional dispersive term appears compared to the previous system. From the LU Boussinesq system, one can also derive a stochastic version of the Kordeveg-De Vries (KdV) equation, as developed in appendix C.

## 3.3 Discussion on conserved quantities

In this section, we discuss the conservation of the following quantities: mass, momentum and (mechanical) energy. We consider a bounded horizontal domain $S_H$ and assume that the noise term $\overline{\sigma \circ d\mathbf{B}_t^H}$ is zero on the boundary $\partial S_H$. In particular, $\overline{\sigma \circ d\mathbf{B}_t^H}$ is periodic.

**Mass:** We regard mass conservation first. Notice that the total mass $m$ fulfils $m = \int_{S_H} \int_0^h \rho dz dS_H$, where $S_H$ is the horizontal domain – which can be 1D or 2D depending on the considered problem – and $h = (1 + \epsilon\eta)$. Assuming that $\rho = \rho_0$ is constant, we get $m = \rho_0 h (|S_H| + \epsilon \int_{S_H} \eta dS_H)$. Hence, $d_t m \propto d_t (\int_{S_H} \eta dS_H)$.

Moreover, in the three models studied in this paper, the evolution equation of the surface elevation $\eta$ remains unchanged and reads

$$d_t \eta + \epsilon (\bar{u}^* \cdot \nabla_H) \eta dt + \Upsilon^{1/2} \epsilon (\overline{\sigma \circ d\mathbf{B}_t^H} \cdot \nabla_H) \eta = -h (\nabla_H \cdot \bar{u}^* dt + \Upsilon^{1/2} \nabla_H \cdot \overline{\sigma \circ d\mathbf{B}_t^H}), \quad (3.68)$$

that is to say, in conservative form,

$$d_t \eta + \nabla_H \cdot (\bar{u}^* h) dt + \Upsilon^{1/2} \nabla_H \cdot (\overline{\sigma \circ d\mathbf{B}_t^H} h) = 0. \quad (3.69)$$

Thus, integrating over the horizontal domain $S_H$ and using the divergence theorem, equation (3.69) yields

$$d_t \left( \int_{S_H} \eta \, dS_H \right) + \int_{\partial S_H} \left( \bar{u}^* h \, dt + \Upsilon^{1/2} \overline{\sigma \circ d\mathbf{B}_t^H} h \right) \cdot dn_{S_H} = 0 \quad (3.70)$$

Therefore, under suitable boundary conditions – which are periodic with a 1D domain $S_H$ in our study – the horizontal integral of the surface elevation $\int_{S_H} \eta \, dS_H$ is conserved, and consequently so is the total mass.

**Momentum:** To investigate momentum conservation, we define the total momentum $\mathbf{p} = \int_{S_H} \int_0^h \rho \bar{u} dz dS_H = \rho_0 \int_{S_H} h \bar{u} \, dS_H$. Then, we derive its evolution equation: starting from the LU Serre-Green-Naghdi model, we
get

$$
\begin{align*}
d_t(h\bar{u}) &= h d_t \bar{u} + \bar{u} d_t h = h d_t \bar{u} + \epsilon \bar{u} d_t \eta \\
&= -\epsilon(h\bar{u}^* \cdot \nabla_H)\bar{u} dt - \Upsilon^{1/2} \epsilon(h\overline{\sigma \circ dB_t^H} \cdot \nabla_H)\bar{u} - h\nabla_H \eta dt + \epsilon\beta^2 \nabla_H \left(\frac{h^3}{3}(d\bar{G})\right) \\
&\quad - \epsilon\bar{u} \left(\nabla_H \cdot (\bar{u}^* h) dt + \Upsilon^{1/2} \nabla_H \cdot (\overline{\sigma \circ dB_t^H} h)\right) \\
&= -\epsilon \nabla_H \cdot (h\bar{u} \otimes \bar{u}^*) dt - \epsilon \Upsilon^{1/2} \nabla_H \cdot (h\bar{u} \otimes \overline{\sigma \circ dB_t^H}) - \underbrace{\frac{1}{2\epsilon} \nabla_H h^2}_{=\nabla_H \left(\eta + \frac{\epsilon\eta^2}{2}\right)} dt + \epsilon\beta^2 \nabla_H \left(\frac{h^3}{3}(d\bar{G})\right). \tag{3.71}
\end{align*}
$$

Similarly as for the mass, for $i \in \{x, y\}$ in 2D, and for $i = x$ in 1D, equation (3.71) yields

$$
d_t \left( \int_{S_H} h \bar{u}_i \, dS_H \right) + \epsilon \int_{\partial S_H} \left( h \bar{u}_i \bar{u}^* \, dt + \Upsilon^{1/2} h \bar{u}_i \overline{\sigma \circ d B_t^H} \right) \cdot \mathbf{d}n_{S_H} + \int_{\partial S_H} \left( \left( \eta + \frac{\epsilon \eta^2}{2} \right) dt - \frac{\epsilon \beta^2 h^3}{3} \overline{dG} \right) \mathbf{d}(\mathbf{n}_{S_H})_i = 0, \quad (3.72)
$$

using the gradient and the divergence theorems. Remind that $d\bar{G}$ is defined as

$$
d\bar{G}(x, t) = \bar{D}_t^H (\nabla_H \cdot \bar{u}) - \epsilon \nabla_H \cdot (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d B_t^H}) \nabla_H \cdot \bar{u}. \quad (3.73)
$$

Again, under suitable boundary conditions – i.e periodic in our work – it is immediate that

$$
\epsilon \int_{\partial S_H} \left( h \bar{u}_i \bar{u}^* dt + \Upsilon^{1/2} h \bar{u}_i \overline{\sigma \circ d \bar{B}_t^H} \right) \cdot \mathbf{d}n_{S_H} + \int_{\partial S_H} \left( \eta + \frac{\epsilon \eta^2}{2} \right) dt \, \mathbf{d}(\mathbf{n}_{S_H})_i = 0. \quad (3.74)
$$

Consequently, one has

$$
d_t \left( \int_{S_H} h \bar{u}_i \, dS_H \right) = \frac{\epsilon \beta^2}{3} \int_{\partial S_H} h^3 d\bar{G} \, d(\mathbf{n}_{S_H})_i. \quad (3.75)
$$

The water height $h$ being periodic by assumption, it is enough to show that $d\bar{G}$ is periodic. In the LU Saint-Venant model, the RHS of equation (3.75) is completely neglected, which is equivalent to assuming $d\bar{G} = 0$. In such case, momentum conservation is immediate. In the LU Boussinesq model, $d\bar{G}$ is approximated as

$$
d\bar{G} = d_t(\nabla_H \cdot \bar{u}) := d\bar{G}_B.
$$

Since $\bar{u}$ is periodic, $\nabla_H \cdot \bar{u}$ also is, as long as the 1st order space derivatives of $\bar{u}$ are well defined. Then $d_t(\nabla_H \cdot \bar{u})$ is periodic as well – as long as this term is well-defined – which proves momentum conservation. Regarding the LU Serre-Green-Naghdi model, we use the “full” equation (3.73) on $d\bar{G}$, namely

$$
\begin{align*}
d\bar{G} = d\bar{G}_{SGN} := \bar{D}_t^H (\nabla_H \cdot \bar{u}) - \epsilon \nabla_H \cdot (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d B_t^H}) \nabla_H \cdot \bar{u} \\
&= d\bar{G}_B + \epsilon (\bar{u}^* \cdot \nabla_H)(\nabla_H \cdot \bar{u}) dt + \Upsilon^{1/2} \epsilon (\overline{\sigma \circ d B_t^H} \cdot \nabla_H)(\nabla_H \cdot \bar{u}) \\
&\quad - \epsilon \nabla_H \cdot (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d B_t^H}) \nabla_H \cdot \bar{u}.
\end{align*}
$$

By similar arguments, the new terms on the RHS are periodic since the 2nd order space derivatives of $\bar{u}$ are,
as long as they are well-defined. Hence, $d\bar{G}_{SGN}$ is periodic as well, that is momentum is conserved.

**Energy:** For shallow water models – in particular, the LU Saint-Venant model – the total energy $E_{SW}$ is
defined as follows (using dimensioned velocities $\bar{u}$ and water height $h$),

$$
\begin{align}
E_{SW} &= \int_{S_H} \int_0^h \frac{1}{2} \rho_0 \|\mathbf{u}\|^2 dz \, dS_H + \int_{S_H} \int_0^h \rho_0 gz \, dz \, dS_h \nonumber \\
&= \frac{\rho_0}{2} \int_{S_H} h \|\bar{u}\|^2 dS_H + \frac{\rho_0 g}{2} \int_{S_H} h^2 dS_H, \tag{3.76}
\end{align}
$$
where the two terms on the RHS respectively correspond to the kinetic and potential energies. Scaling $\bar{u}$ and $h$ as before, we find the equation on the following rescaled energy equation [32, 33]

$$E_{SW} = \frac{\epsilon^2}{2} \int_{S_H} h \|\bar{u}\|^2 dS_H + \frac{1}{2} \int_{S_H} h^2 dS_H = \frac{\epsilon^2}{2} (h\bar{u}, \bar{u})_{L^2(S_H, \mathbb{R}^2)} + \frac{1}{2} \|h\|_{L^2(S_H, \mathbb{R})}^2 := E_c + E_p, \quad (3.77)$$

denoting $E_c$ and $E_p$ the scaled total kinetic and potential energies. We also denote $e_c = \frac{\epsilon^2}{2} h \|\bar{u}\|^2$ and $e_p = \frac{1}{2} h^2$. Now, we derive the evolution equation of this energy in the LU Saint-Venant model, using Einstein's notation on $i$,

$$\begin{aligned} d_t e_c &= \frac{\epsilon^2}{2} \bar{u} \cdot d_t(h\bar{u}) + \frac{\epsilon^2}{2} h\bar{u} \cdot d_t \bar{u} \\ &= -\frac{\epsilon^3}{2} \sum_{j \in J} \bar{u}_i \partial_j (h\bar{u}_i \bar{u}_j^*) dt - \frac{\epsilon^3 \Upsilon^{1/2}}{2} \sum_{j \in J} \bar{u}_i \partial_j (h\bar{u}_i (\overline{\sigma \circ d B_t^H})_j) - \frac{\epsilon}{2} h\bar{u}_i \partial_i h dt \\ &\quad - \frac{\epsilon^3}{2} \sum_{j \in J} h\bar{u}_i \bar{u}_j^* \partial_j \bar{u}_i dt - \frac{\epsilon^3 \Upsilon^{1/2}}{2} \sum_{j \in J} h\bar{u}_i (\overline{\sigma \circ d B_t^H})_j \partial_j \bar{u}_i - \frac{\epsilon}{2} h\bar{u}_i \partial_i h dt \\ &= -\frac{\epsilon^3}{2} \nabla_H \cdot (h\bar{u}_i^2 \bar{u}^*) dt - \frac{\epsilon^3}{2} \Upsilon^{1/2} \nabla_H \cdot (h\bar{u}_i^2 (\overline{\sigma \circ d B_t^H})_j) - \epsilon h\bar{u}_i \partial_i h dt, \end{aligned}$$

and

$$d_t e_p = h \cdot d_t h = -\epsilon h \partial_i (h \bar{u}_i^*) dt - \epsilon \Upsilon^{1/2} h \partial_i (h (\overline{\sigma \circ d B_t^H})_i),$$

where $J = \{x\}$ is the problem is 1D and $J = \{x, y\}$ if it is 2D. Thus, the quantity $e = e_c + e_p$ fulfils

$$\begin{aligned} \nabla \cdot {}_t e = & -\frac{\epsilon^3}{2} \nabla_H \cdot (h \|\bar{u}\|^2 \bar{u}^*) dt - \frac{\epsilon^3}{2} \Upsilon^{1/2} \nabla_H \cdot (h \|\bar{u}\|^2 \overline{\sigma \circ d B_t^H}) - \epsilon \nabla_H \cdot (h^2 \bar{u}) dt \\ & + \frac{\epsilon \Upsilon}{4} h \nabla_H \cdot (h \bar{u}_s) dt - \epsilon \Upsilon^{1/2} h \nabla_H \cdot (h \overline{\sigma \circ d B_t^H}). \end{aligned} \quad (3.78)$$

Integrating over the domain $S_H$, using the divergence theorem and periodic boundary conditions, we get by integration by parts the evolution equation of the total Saint-Venant energy $E_{SV}$,

$$\begin{aligned} d_t E_{SV} &= \int_{S_H} \left[ -\frac{\epsilon \Upsilon}{8} \bar{u}_s \cdot \nabla_H h^2 dt + \frac{\epsilon}{2} \Upsilon^{1/2} \overline{\sigma \circ d B_t^H} \cdot \nabla_H h^2 \right] dS_H \\ &= \int_{S_H} \left[ \frac{\epsilon \Upsilon}{8} h^2 \nabla_H \cdot \bar{u}_s dt - \frac{\epsilon}{2} \Upsilon^{1/2} h^2 \nabla_H \cdot \overline{\sigma \circ d B_t^H} \right] dS_H. \end{aligned} \quad (3.79)$$

Consequently, for the LU Saint-Venant model – that is assuming $d\bar{G} = 0$ – energy conservation is enforced by choosing the noise term $\sigma \circ d B_t$ such that $\nabla_H \cdot \sigma \circ d B_t^H = \nabla_H \cdot \bar{u}_s = 0$. We denote this assumption (DF-BHNISD), standing for “divergence free barotropic horizontal noise and Itô-Stokes drift”. However, in 1D problems, this condition does not make much physical sense since it is equivalent to considering a constant horizontal noise over space. For this reason, energy conservation is not ensured in our numerical simulations, since they were performed with more general noises which do not fulfil this assumption. Nevertheless, we anticipate that performing 2D test simulations of this stochastic Saint-Venant equation with non trivial divergence free noise would lead to numerical energy conservation results.

For the LU Serre-Green-Naghdi model, the energy is rather defined as [32, 34]

$$E_{SGN} = \frac{\epsilon^2}{2} (h\bar{u}, \bar{u})_{L^2(S_H, \mathbb{R}^2)} + \frac{1}{2} \|h\|_{L^2(S_H, \mathbb{R})}^2 + \frac{\epsilon^3 \beta^2}{6} (h^3 \nabla_H \cdot \bar{u}, \nabla_H \cdot \bar{u})_{L^2(S_H, \mathbb{R})}. \quad (3.80)$$

Using the same notations $e_c$ and $e_p$ as before, and defining $e_{pv} = \frac{1}{6} h^3 (\nabla_H \cdot \bar{u})^2$ and $e = e_c + e_p + e_{pv}$, one has
similarly

$$d_t E_{SGN} = \int_{S_H} \left[ -\frac{\epsilon^3 \beta^2}{3} h^3(\nabla_H \cdot \bar{u}) d\bar{G} + \frac{\epsilon^3 \beta^2}{2} h^2 d_t h(\nabla_H \cdot \bar{u})^2 + \frac{\epsilon^3 \beta^2}{3} h^3(\nabla_H \cdot \bar{u}) d_t(\nabla_H \cdot \bar{u}) \right] dS_H, \quad (3.81)$$

using that $\nabla_H \cdot \overline{\sigma \circ d\bar{B}_t^H} = \nabla_H \cdot \bar{u}_s = 0$. Now, computing the first term in the integrand yields

$$-\frac{\epsilon^3 \beta^2}{3} h^3 (\nabla_H \cdot \bar{u}) d\bar{G} = -\frac{\epsilon^3 \beta^2}{3} h^3 (\nabla_H \cdot \bar{u}) d_t (\nabla_H \cdot \bar{u}) - \frac{\epsilon^4 \beta^2}{3} h^3 (\nabla_H \cdot \bar{u}) (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) \cdot \nabla_H (\nabla_H \cdot \bar{u}) \\ + \frac{\epsilon^4 \beta^2}{3} h^3 (\nabla_H \cdot \bar{u})^2 \nabla_H \cdot (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}), \quad (3.82)$$

that is

$$-\frac{\epsilon^3 \beta^2}{3} h^3 (\nabla_H \cdot \bar{u}) d\bar{G} + \frac{\epsilon^3 \beta^2}{3} h^3 (\nabla_H \cdot \bar{u}) d_t (\nabla_H \cdot \bar{u}) = -\frac{\epsilon^4 \beta^2}{6} h^3 (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) \cdot \nabla_H (\nabla_H \cdot \bar{u})^2 \\ + \frac{\epsilon^4 \beta^2}{3} h^3 (\nabla_H \cdot \bar{u})^2 \nabla_H \cdot (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}). \quad (3.83)$$

In addition, the second term in the integrand is

$$\frac{\epsilon^3 \beta^2}{2} h^2 d_t h (\nabla_H \cdot \bar{u})^2 = -\frac{\epsilon^4 \beta^2}{2} h^2 (\nabla_H \cdot \bar{u})^2 \nabla_H \cdot (h \bar{u}^* dt + \Upsilon^{1/2} h \overline{\sigma \circ d\bar{B}_t^H}) \quad (3.84) \\ = -\frac{\epsilon^4 \beta^2}{2} h^3 (\nabla_H \cdot \bar{u})^2 \nabla_H \cdot (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) - \frac{\epsilon^4 \beta^2}{6} (\nabla_H \cdot \bar{u})^2 (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) \cdot \nabla_H h^3.$$

Consequently,

$$d_t E_{SGN} = \frac{\epsilon^4 \beta^2}{6} \int_{S_H} \left[ h^3 (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) \cdot \nabla_H (\nabla_H \cdot \bar{u})^2 \right] dS_H \\ + \frac{\epsilon^4 \beta^2}{6} \int_{S_H} \left[ h^3 (\nabla_H \cdot \bar{u})^2 \nabla_H \cdot (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) \right] dS_H \\ + \frac{\epsilon^4 \beta^2}{6} \int_{S_H} \left[ (\nabla_H \cdot \bar{u})^2 (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) \cdot \nabla_H h^3 \right] dS_H \\ = \frac{\epsilon^4 \beta^2}{6} \int_{S_H} \nabla_H \cdot \left[ h^3 (\nabla_H \cdot \bar{u})^2 (\bar{u}^* dt + \Upsilon^{1/2} \overline{\sigma \circ d\bar{B}_t^H}) \right] dS_H = 0, \quad (3.85)$$

using the divergence theorem and the periodic boundary conditions again. Notice that no assumptions where made in addition to the one for the Saint-Venant model, that is (DF-BHNISD). As before, this assumption leads to a space constant noise in the 1D case, therefore it is anticipated that the energy is not conserved in our simulations. Moreover, the previous calculations show that (DF-BHNISD) is not enough to enforce energy conservation in the LU Boussinesq model for both the energies $E_{SW}$ and $E_{SGN}$, which is coherent with the deterministic Boussinesq model.

# 4 Numerical simulations

In this section, we present some numerical simulations we made to test the three models derived. The Julia code that we produced is based on the work of Vincent Duchêne and Pierre Navarro, who proposed a variety of wave models implementations in the deterministic setting – see the documentation in the following link: http://waterwavesmodels.github.io/WaterWaves1D.jl/dev/ [35]. These models are essentially based on pseudo-spectral resolution methods, which justifies the use of periodic boundary conditions. We adapted their numerical framework to the stochastic case, introducing implementations of the noise terms and the ISDs for this purpose.

Regarding the purely stochastic aspects, we consider noises with wave spatial structure. This is justified by the shape of solutions found in [3]: considering a constant noise a first, the authors showed that the system admits
progressive wave solutions. Then, they extend the analysis to systems where the noise is itself a progressive wave. In the end, the 1D noise we consider is the following,

$$ \sigma(x)dW_t = A \cos(kx)d\beta_t^1 + A \sin(kx)d\beta_t^2, \qquad (4.1) $$

where A denotes the amplitude of the noise and $d\beta_t^1, d\beta_t^2$ denote Brownian motions. Using a Box-Muller argument, this shape is equivalent to

$$ \sigma(x)dW_t = A \cos(kx + \phi_t)d\beta_t, \qquad (4.2) $$

where $\phi_t$ is a uniformly distributed random phase on $(-\pi, \pi)$, such that for all $t \neq s$, $\phi_t$ and $\phi_s$ are independent. Additionally, $d\beta_t$ is a Brownian motion. In our simulations, the noise wave number is set to $k = 2\pi/100$, and the noise amplitude may take the following values: $A = 0.001$, $A = 0.005$ or $A = 0.01$. Notice that the dimensioned noise scales like $A\epsilon\sqrt{gh_0}$ as a consequence. The value of wave number $k$ is chosen to be at least one order of magnitude smaller than the typical wave number of the deterministic wave. This is because our simulations showed that noise terms with too small space scale oscillations lead to numerical instability, and enables us to further discuss the presence of an additive noise term in the water elevation dynamics. The values of amplitude where chosen to be much smaller than the typical height of the wave. Numerically, we have observed that $A = 0.001$ yields slight perturbations of the deterministic waves – that is typical realisations of each LU model is similar to its deterministic counterpart – while $A = 0.01$ induces a more “noisy” dynamics – that is typical realisations are essentially noise driven. In addition, we chose $A = 0.005$ as an enlightening intermediate case.

Our tests are based on computing the evolution of the deformation surface $\eta$, with a “heap of water” type initial condition. Namely, the initial surface deformation is set to be $\eta(x, t = 0) = \exp(-x^4)$, while the initial velocity is set to $u(x, t = 0) = 0$. All of our tests were performed on a numerical 1D tank $[-L, L]$ with $L = 50$, which is discretised with $N = 2^{11}$ spatial points. The timestep is chosen to be $dt = 0.005s$, and we assume periodic boundary conditions. To enforce these conditions on the noise terms as well, we multiply them the function $s_\alpha(x) = \exp\left(\frac{1}{\alpha^2}\left(1 - \frac{1}{1-(x/L)^2}\right)\right)$, with $\alpha = 10$, in order to make them vanish on the boundary. The initial condition on $\eta$ and the profile of $s_\alpha$ are given on figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3: Initial surface deformation (left) and boundary conditions enforcement function $s_\alpha$ (right, $\alpha = 10$).</figcaption></figure>

Moreover, we propose two sets of parameters for testing our models

* $\beta = 0.01$ and $\epsilon = 0.1$ ($\mathcal{P}_1$), so that the scaling conditions associated with the Shallow Water model hold (and in particular, they hold for the Boussinesq model as well). Therefore, the three models should give qualitatively similar results.

* $\beta = 0.1$ and $\epsilon = 0.1$ ($\mathcal{P}_2$), so that the scaling conditions of the Shallow Water and the Boussinesq models do not hold. Thus, these models should give qualitatively distinct – yet not dramatically different – results.
Additionally, we will compare our stochastic models to the deterministic one dimensional Water Waves model.

Furthermore, two options are available for simulating our models, depending on whether we consider a “true” Stratonovitch noise or if we rewrite it in $It\bar{o}$ form. In the former case, we may adopt a stochastic Euler-Heun approach as developed in [36], where the computation of the stochastic diffusion becomes implicit. In the latter case however, one needs to give an analytical expression for the stochastic diffusion. Even though one can compute them in the LU Saint-Venant and the LU Boussinesq systems, this correction term becomes extremely complex in the LU Serre-Green-Naghdi model due to the term $\overline{dG}$. Therefore, we rather adopted the first approach using the stochastic Euler-Heun method.

Moreover, for numerical stability reasons, we used a stochastic version of the order 4 Runge-Kutta algorithm (RK4), rather than the (simpler) Euler-Maruyama algorithm. For instance, the (deterministic) Saint-Venant equations are known to be dramatically more stable when solved with the order 4 Runge-Kutta rather than the Euler method. In summary, the solving algorithm we used is essentially the following: we treat the bounded variation term as in the classical RK4 method, and the martingale term as in the stochastic Euler-Heun method. Such approach has been studied in more details in [37, 38, 39].

Furthermore, since water elevation equations on $\eta$ are the same in the three models we study – regarding for example equation (3.68) – one may notice the presence of the term

$$ \Upsilon^{1/2}h\nabla_H \cdot (\overline{\sigma \circ d\mathbf{B}_t^H}) = \Upsilon^{1/2}\nabla_H \cdot (\overline{\sigma \circ d\mathbf{B}_t^H}) + \epsilon\Upsilon^{1/2}\eta\nabla_H \cdot (\overline{\sigma \circ d\mathbf{B}_t^H}). $$

This shows the existence of an additive noise term $\Upsilon^{1/2}\nabla_H \cdot (\overline{\sigma \circ d\mathbf{B}_t^H})$ in the water elevation dynamics. Its effect is illustrated by Figure 4: due to this term, the surface elevation is not flat “away from the wave”, which is a strong difference compared to the deterministic setting. To facilitate the comparison between our models and their deterministic counterparts, we chose to disregard this additive noise term in our simulations. This can be interpreted as a filtering of the lower wave numbers – i.e. large scale dynamics.

<figure><img src="image_5.png" /><figcaption>Figure 4: Realisations of the LU Saint-Venant model, with additive noise (left) and without (right). We only plot the solution over the domain $[-20, 20]$, at $t = 5s$. Parameter set: ($\mathcal{P}_1$) – Wave number: $k = 2\pi/100$ – Amplitude: 0.001.</figcaption></figure>

## 4.1 Qualitative analysis on the effect of the noise

In this subsection, we give some qualitative insights about the effect of the noise on the dynamics of the system. For this purpose, we compare each deterministic solution to a solution of the associated LU interpretation, disregarding the effect of the additive noise term previously mentioned.
### 4.1.1 Parameters ($\mathcal{P}_1$)

Considering the parameters ($\mathcal{P}_1$), we compared the deterministic Saint-Venant, Boussinesq and Serre-Green-Naghdi models (blue curves) to realisations of their associated LU models (orange curves), using the 1D Water Waves as a reference (green curves). The noise we chose has the shape of a stationary wave, with wave number $k = 2\pi/10$. Figure 5 shows realisations of these models at $t = 5s$, for different noise amplitudes.

<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison between the surface deformation of the deterministic Saint-Venant (1st row), Boussinesq (2nd row) and Serre-Green-Naghdi (3rd row) models and realisations of their LU interpretations. Parameter set: ($\mathcal{P}_1$) – Wave number: $k = 2\pi/100$ – Amplitude, from left to right: 0.001, 0.005 and 0.01.</figcaption></figure>

We observe that deterministic solutions are close to the Water Waves solution. This is expected regarding the scaling of $\beta = 0.01 \ll 1$. Moreover, the LU systems seem to converge to their respectively associated deterministic solution for a vanishing noise, thus giving consistency to the LU interpretation of the wave equations. In both cases, the noise tends to break the spatial symmetry of the wave. This stems from the shape of the noise, which is a sum of a symmetric function – $A \cos(kx)d\beta_t^1$ – and an antisymmetric function – $A \sin(kx)d\beta_t^2$.

### 4.1.2 Parameters ($\mathcal{P}_2$)

The symmetry breaking mentioned in the previous subsection is also observed in 6, where we use parameters ($\mathcal{P}_2$). However, in this case, the parameter $\beta = 0.1$ does not match the validity conditions of the Saint-Venant equations since $\beta \sim \epsilon$. Therefore, the Water Waves solution is expected to differ from the Saint-Venant ones, in both deterministic and LU forms. This is observed in figure 6, indeed. In addition, the LU Boussinesq and LU Serre-Green-Naghdi appear to give an interesting variability to their deterministic versions. However, the LU Serre-Green-Naghdi model yields numerical instabilities when choosing larger values of $\beta$ – e.g. $\beta = 1$. This is due to a periodicity default caused by the symmetry breaking of the wave, which does not exist in the deterministic setting. To tackle this issue, we think of investigating numerical models in conservative form – that is computing the water elevation/momentum variables $(\eta, h\bar{u})$ rather than the water elevation/momentum variables $(\eta, \bar{u})$. For similar reasons, our implementation does not enjoy conservation of physical quantities
discussed earlier. Although it is relatively simple to translate our LU Saint-Venant algorithm into conservative form – which then conserves mass, momentum and energy up to machine accuracy – our LU Boussinesq and LU Serre-Green-Naghdi algorithms are more challenging to adapt due to the presence of the term $\overline{dG}$. We expect such numerical method to be more stable than the currently used one, and to allow investigating how multi-scale location uncertainty affects the waves dynamics. This is subject to further work.

<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison between the surface deformation of the deterministic Saint-Venant (1st row), Boussinesq (2nd row) and Serre-Green-Naghdi (3rd row) models and their LU interpretations. Parameter set: ($\mathcal{P}_2$) – Wave number: $k = 2\pi/100$ – Amplitude, from left to right: 0.001, 0.005 and 0.01. The bottom right figure illustrates that the LU Serre-Green-Naghdi model is unstable for large enough noise amplitudes.</figcaption></figure>

## 4.2 Numerical estimation of the noise-induced spreading

In this section, we analyse the first and second order statistics of each LU wave model in the setting described above, at time $t = 5s$. Again, we will study these models in the sets of parameters ($\mathcal{P}_1$) for the LU Saint-Venant model, and ($\mathcal{P}_2$) for the LU Boussinesq and Serre-Green-Naghdi models. The LU models statistics we analyse are computed with 130 realisations of each stochastic model.

A spreading arises from LU Saint-Venant model, which appears to grow linearly with the amplitude of the noise – see Fig 7. Also, is concentrated in the “upstream” of the wave, which is expected since it would be physically irrelevant for the wave to affect a region where it has not passed. In addition, the means and the standard deviations appear to be space symmetric – up to statistical error – which suggests the distribution is also space symmetric. Moreover, there exists a significant spreading at the peak of the wave for strong enough noise amplitudes.

Regarding the LU Boussinesq and Serre-Green-Naghdi models, the same remarks on the peak spreading and the symmetry of the means and standard deviations apply – see Fig 8. Again, the spreading is concentrated in the “upstream” of the wave, specifically at its peaks and troughs. In particular, the maximum height value of the LU models varies depending on the stochastic realisation. Such observation gives lines of approach for building stochastic models by selecting the noise $\sigma$ to data, using for instance calibration or data assimilation
<figure><img src="image_8.png" /><figcaption>Figure 7: First and second order statistics of the LU Saint-Venant model at $t = 5s$, compared to the associated deterministic model and the water waves one. In addition, an evaluation of the spreading – defined here as 3 times the empirical standard deviation – is given for different values of noise amplitude (orange area). We only plot the solution over the domain $[-20, 20]$. Parameter set: ($P_1$) – Wave number: $k = 2\pi/100$ – Amplitude, from left to right: 0.001, 0.005 and 0.01.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: First and second order statistics of the LU Boussinesq (1st row) and LU Serre-Green-Naghdi (2nd row) models at $t = 5s$, compared to the associated deterministic models and the water waves one. In addition, an evaluation of the spreading – defined here as 3 times the empirical standard deviation – is given for different values of noise amplitude (orange area). We only plot the solution over the domain $[-20, 20]$. Parameter set: ($P_2$) – Wave number: $k = 2\pi/100$ – Amplitude, from left to right: 0.001, 0.005 and 0.01.</figcaption></figure>

techniques.

Although the Serre-Green-Naghdi and the Boussinesq models give similar results for the parameters ($P_2$), we expect that much more differences would be observed between the two models with $\beta = 1$ and $\epsilon = 0.1$. As mentioned before, wave models written in conservative form should be more stable and allow to perform such tests in this configuration.

# 5 Conclusion and discussion

In this work, we investigated the stochastic representation of several shallow water coastal wave models within the LU framework. These stochastic models maintain the same conservation properties formally and thus exhibit physical consistency with their deterministic counterparts. We demonstrated numerically that they induce a pathwise symmetry breaking, accompanied by a restoration of this symmetry in law, which is a property that should be expected in representing turbulent effects. For low noise amplitudes, the studied models have shown to converge toward the deterministic solutions. For large amplitudes, we need to transition to a numerical scheme in conservation form. This will be the subject of a future study.
# A Appendix 1: Leibniz formula

We give in the following the expression of the Leibniz integral rule for a stochastic process. We want to evaluate the derivative of an integral of the form

$$ \int_{a(x,t)}^{b(x,t)} f(x, z, t) dz, $$

with respect to the space, $x$, and time, $t$, and where $a(x, z, t)$, $b(x, z, t)$ and $f(x, z, t)$ are continuous $C^2$-semimartingale. For the space variable, no time derivative is at play and the usual formula stands:

$$ \frac{\partial}{\partial x} \int_{a(x,t)}^{b(x,t)} f(x, z, t) dz = \int_{a(x,t)}^{b(x,t)} \partial_x f(x, z, t) dz + \partial_x b f(x, b(x,t), t) - \partial_x a f(x, a(x,t), t). $$

Differentiation with respect to time, involves now a stochastic integration. Understanding the stochastic integral in the Stratonovich setting, we have,

$$ d_t \int_{a(x,t)}^{b(x,t)} f(x, z, t) dz = \int_{a(x,t)}^{b(x,t)} d_t f(x, z, t) dz + d_t b f(x, b(x,t), t) - d_t a f(x, a(x,t), t). $$

The key argument of proof consists in defining the functions $F(x, y, t) = \int_{z_1}^y f(x, z, t) dz$ and $G(x, a, b, t) = \int_{a(x,t)}^{b(x,t)} f(x, z, t) dz$, then link them with a functional relation, that is

$$ \begin{aligned} G(x, a, b, t) &= \int_{0}^{b(x,t)} f(x, z, t) dz - \int_{0}^{a(x,t)} f(x, z, t) dz \\ &= F(x, b(x,t), t) - F(x, a(x,t), t). \end{aligned} $$

Now we introduce the following chain rule for Stratonovich calculus,

**Theorem A.1 (Generalized Itō's formula – Stratonovich form).** Let $\theta(\boldsymbol{x}, t)$, $\boldsymbol{x} \in \Omega$ be a continuous $C^3$-process and a continuous $C^2$-semimartingale, let $\boldsymbol{X}_t$ be a continuous semimartingale with values in $\Omega$. Then, the following formula is satisfied :

$$ d\theta(\boldsymbol{X}_t, t) = d_t\theta(\boldsymbol{X}_t, t) + \frac{\partial\theta}{\partial x_i}(\boldsymbol{X}_t, t) \circ dX_t^i. \quad (A.1) $$

Upon applying this on $G(x, a, b, t)$, we obtain

$$ \begin{aligned} d_t G &= d_t F(x, b(x,t), t) + \partial_b F(x, b(x,t), t) d_t b - d_t F(x, a(x,t), t) - \partial_b F(x, b(x,t), t) d_t a \\ &= \int_{a(x,t)}^{b(x,t)} d_t f(x, z, t) dz + d_t b f(x, b(x,t), t) - d_t a f(x, b(x,t), t), \end{aligned} $$

To interpret the stochastic integral in the Itō setting, we would need to adapt the previous chain rule as follows,

**Theorem A.2 (Generalized Itō's formula – Itō form).** Let $\theta(\boldsymbol{x}, t)$, $\boldsymbol{x} \in \Omega$ be a continuous $C^2$-process and a continuous $C^1$-semimartingale, let $\boldsymbol{X}_t$ be a continuous semimartingale with values in $\Omega$. Then, $\theta(\boldsymbol{X}_t, t)$ is a continuous semimartingale satisfying

$$ d\theta(\boldsymbol{X}_t, t) = d_t\theta(\boldsymbol{X}_t, t) + \frac{\partial\theta}{\partial x_i}(\boldsymbol{X}_t, t)dX_t^i + \frac{1}{2}\frac{\partial^2\theta}{\partial x_i\partial x_j}(\boldsymbol{X}_t, t)d\langle X^i, X^j \rangle_t + d\langle \frac{\partial\theta}{\partial x_i}(\boldsymbol{X}, \cdot), X^i \rangle_t. \quad (A.2) $$
Applying it to the process $G(x, a, b, t)$, we obtain

$$d_t G = \int_{a(x,t)}^{b(x,t)} d_t f(x, z, t) dz + d_t b f(x, b(x, t), t) - d_t a f(x, b(x, t), t) + \\ d\langle f(x, b(x, t), t), b\rangle_t dt - d\langle f(x, a(x, t), t), a\rangle_t dt + \\ \frac{1}{2} d\langle b, b\rangle_t \partial_{b^2}^2 F(x, b(x, t), t) - \frac{1}{2} d\langle a, a\rangle_t \partial_{a^2}^2 F(x, a(x, t), t) dt. \quad (A.3)$$

This latter includes several additional quadratic variation terms, and is more cumbersome to use in formal developments. However, this expression is necessary to access to the mathematical expectation.

## B Appendix 2: Quadratic covariation

In stochastic calculus, the quadratic covariation (or cross-variance) of two processes X and Y is defined as

$$\langle X, Y \rangle_t = \lim_{n \to 0} \sum_{i=1}^{p_n} (X_i^n - X_{i-1}^n)(Y_i^n - Y_{i-1}^n), \quad (B.1)$$

where $0 = t_0^n < t_1^n < \dots < t_{p_n}^n = t$ is a partition of the interval $[0, t]$. The previous limit, if it exists, is defined in the sense of convergence in probability.

Let X and Y two continuous semimartingales, defined as $X_t = X_0 + A_t + M_t$, $Y_t = Y_0 + B_t + N_t$, $M, N$ being martingales and A, B finite variation processes. In this context, their quadratic covariation (B.1) exists, and is given by

$$\langle X, Y \rangle_t = \langle M, N \rangle_t. \quad (B.2)$$

In particular, the quadratic variation of a standard Brownian motion B (as a martingale) is $\langle B \rangle_t := \langle B, B \rangle_t = t$, by definition.

Quadratic covariances play an important role in stochastic calculus, as they arise in Itō's lemma, which can be interpreted as a stochastic chain rule. In particular, these terms are involved in the Itō integration by parts formula,

$$d_t(XY) = Xd_tY + Yd_tX + d\langle X, Y \rangle_t, \quad (B.3)$$

and in Itō's isometry, expressing the covariance of two Itō integrals: let f and g two predictable processes such that $\int_0^t f^2 d\langle M, M \rangle_s$ and $\int_0^t g^2 d\langle N, N \rangle_s$ are finite, then

$$\mathbb{E}[(\int_{0}^{t} f dM_{s})(\int_{0}^{t} g dN_{s})] = \mathbb{E}[\int_{0}^{t} f g d\langle M, N \rangle_{s}]. \quad (B.4)$$

## C Appendix 3: LU Kordeveg - De Vries equation

To derive the LU interpretation of the Kordeveg - De Vries equation (KdV), we will adapt a standard procedure to the stochastic case [40]. For simplification purpose, we will assume homogeneity in the transverse direction and consider in the following a 1D version of the Boussinesq model. Assuming that the bottom is flat, the 1D Boussinesq model reads

$$D_t^H \eta = -h \left( \partial_x (\bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s) dt + \Upsilon^{1/2} \partial_x \overline{\sigma \circ d \mathbf{B}_t^x} \right), \quad (C.1a)$$

$$d_t \bar{u} + \epsilon \left( \bar{u} - \frac{1}{2} \Upsilon \epsilon \bar{u}_s \right) \partial_x \bar{u} dt + \Upsilon^{1/2} \epsilon \overline{\sigma \circ d \mathbf{B}_t^x} \partial_x \bar{u} + \partial_x \eta dt - \epsilon \beta^2 \left( \frac{h^2}{3} \partial_{xx}^2 d_t \bar{u} \right) = \mathcal{O}(\beta^4, \epsilon \beta^4) \quad (C.1b)$$
The solutions we consider for the above Boussinesq system are assumed to be waves, and we apply the following change of variable: $\xi = x - t + \varphi_t$, where $\varphi_t = f(B_t)$ is a random phase that does not depends on $x$. The wave shape is also assumed to change on a large temporal scale $\tau = \epsilon t$. The depth averaged horizontal velocity $u(\xi, \tau, \varphi)$ and the surface elevation $\eta(\xi, \tau, \varphi)$ are assumed to be smooth functions. The noise $\sigma_0 dB_t^x(\xi, \tau)$ is assumed to be homogeneous, and thus is associated to a constant variance tensor and a zero ISD. In this new formalism the surface elevation equation reads

$$-\partial_{\xi}\eta dt + \epsilon\partial_{\tau}\eta dt + \partial_{\varphi}\eta f'(B_t)dB_t + \epsilon u\partial_{\xi}\eta dt + \Upsilon^{1/2}\epsilon \sigma_0 dB_t^x \partial_{\xi}\eta + h(\partial_{\xi}u dt + \Upsilon^{1/2}\partial_{\xi}\sigma_0 dB_t^x) = 0. \quad (C.2)$$

After converting this equation in $It\bar{o}$ form, terms of finite variations (ie “$dt$” terms) and martingale terms (i.e. “$dB_t$” terms) can be rigorously separated by the Biechteller-Delacherie theorem. For the martingale terms we have

$$\partial_{\varphi}\eta f'(B_t)dB_t = -\Upsilon^{1/2}\partial_{\xi}(h\sigma dB_t^x), \quad (C.3)$$

and for the finite variation terms

$$-\partial_{\xi}\eta + \epsilon\partial_{\tau}\eta + \epsilon u\partial_{\xi}\eta - \frac{1}{2}\Upsilon\epsilon^2 a^H \partial_{\xi\xi}\eta + h\partial_{\xi}u = 0. \quad (C.4)$$

Moreover, for the velocity equation we have

$$-\partial_{\xi}u dt + \epsilon\partial_{\tau}u dt + \partial_{\varphi}uf'(B_t)dB_t + \epsilon u\partial_{\xi}u dt + \Upsilon^{1/2}\epsilon \sigma_0 dB_t^x \partial_{\xi}u \\ -\frac{1}{2}\Upsilon\epsilon^2 \bar{a}^{hh}\partial_{\xi\xi}^2 u dt + \partial_{\xi}\eta dt - \frac{1}{h}\epsilon\beta^2\left(\frac{h^3}{3}\partial_{\xi\xi}^2(-\partial_{\xi}u + \epsilon\partial_{\tau}u)\right) dt = \mathcal{O}(\beta^4, \epsilon\beta^4). \quad (C.5)$$

Then, the martingale terms yield

$$\partial_{\varphi}uf'(B_t)dB_t = -\Upsilon^{1/2}\epsilon \sigma_0 dB_t^x \partial_{\xi}u \quad (C.6)$$

and for the finite variation terms

$$-\partial_{\xi}u + \epsilon\partial_{\tau}u + \epsilon u\partial_{\xi}u - \frac{1}{2}\Upsilon\epsilon^2 a^H \partial_{\xi\xi}u + \partial_{\xi}\eta + \frac{1}{h}\epsilon\beta^2\left(\frac{h^3}{3}(\partial_{\xi\xi\xi}^3 u - \epsilon\partial_{\xi\xi\tau}u)\right) = 0. \quad (C.7)$$

Expanding $u$ and $\eta$ in terms of the small parameter $\epsilon$ as $q = q_0 + \epsilon q_1 + \dots$, and identifying the equations term by term of corresponding order, we get for the zero order terms,

$$\partial_{\xi}\eta_0 = \partial_{\xi}u_0 \quad \text{and hence} \quad \eta_0 = u_0.$$

At order $\epsilon$, for a noise of magnitude up to $\Upsilon \sim \mathcal{O}(1)$, we obtain the system

$$-\partial_{\xi}u_1 + \partial_{\tau}\eta_0 + \eta_0\partial_{\xi}\eta_0 + \partial_{\xi}\eta_1 + \frac{1}{3}\partial_{\xi\xi\xi}^3\eta_0 = 0, \\ -\partial_{\xi}\eta_1 + \partial_{\tau}\eta_0 + 2\eta_0\partial_{\xi}\eta_0 + \partial_{\xi}u_1 = 0,$$

from which one gets immediately the classical KdV equation with random phase

$$\partial_{\tau}\eta_0 + \frac{3}{2}\eta_0\partial_{\xi}\eta_0 + \frac{1}{6}\partial_{\xi\xi\xi}^3\eta_0 = 0. \quad (C.8)$$

This equation has the structure of a Burger equation with an additional dispersive term. A modified KdV equation is obtained by considering a stronger noise of amplitude $\Upsilon \sim 1/\epsilon$. In that case the second order terms must be kept in the Taylor development of terms at $\epsilon$ order. Eventually, one obtains the dissipative KdV
equation (with random phase)

$$ \partial_{\tau} \eta_0 + \frac{3}{2} \eta_0 \partial_{\xi} \eta_0 - \frac{1}{2} a^H \partial_{\xi \xi} \eta_0 + \frac{1}{6} \partial_{\xi \xi \xi}^3 \eta_0 = 0. \quad (C.9) $$

The random phase of both KdV equation is determined by (C.3) and (C.6). From the latter, we notice immediately that $f'(B)$ and $\epsilon$ must share the same order for the solution not to be trivial. At order $\epsilon$ from (C.6), one has

$$ \partial_{\varphi} \eta_0 f'(B_t) dB_t = -\Upsilon^{1/2} \partial_{\xi} \eta_0 \sigma d\mathbf{B}_t^{0,x}. $$

A simple way to choose $f$ is to impose $f(B_t) = -k\Upsilon^{1/2}\epsilon k^\sigma B_t$ and $\sigma = k^\sigma$, where $k^\sigma$ is associated to the waves form

$$ \eta_0 = h_0 e^{i(kx - \omega t - \epsilon \Upsilon^{1/2} k k^\sigma B_t)}. \quad (C.10) $$

This solution corresponds to stochastic linear waves solutions as derived in [4]. Note that the derivations above have been done using the strong hypothesis of $u$ and $\eta$ being smooth functions of time – more precisely of finite variation. Looking now for more general stochastic solutions, one considers the following variables

$$ u(\xi, \tau), \quad \sigma \circ d\mathbf{B}_t^H(\xi, \tau), \quad \eta(\xi, \tau), \quad (C.11) $$

which are not differentiable with respect to $\tau$ (i.e. they are semi-martingale stochastic processes). Making the same assumptions and deriving the equations in Stratonovich form in the same way as previously, we obtain two coupled SPDE's,

$$ -\partial_{\xi} \eta dt + \epsilon d_{\tau} \eta + \epsilon u \partial_{\xi} \eta dt + \Upsilon^{1/2} \epsilon \sigma \circ d\mathbf{B}_t \partial_{\xi} \eta + h(\partial_{\xi} u dt + \Upsilon^{1/2} \partial_{\xi} \sigma \circ d\mathbf{B}_t^x) = 0, \quad (C.12) $$

and

$$ -\partial_{\xi} u dt + \epsilon d_{\tau} u + \epsilon u \partial_{\xi} u dt + \Upsilon^{1/2} \epsilon \sigma \circ d\mathbf{B}_t^x \partial_{\xi} u + \partial_{\xi} \eta dt - \frac{1}{h} \epsilon \beta^2 \left( \frac{h^3}{3} \partial_{\xi \xi}^2 (-\partial_{\xi} u + \epsilon \partial_{\tau} u) \right) dt = \mathcal{O}(\beta^4, \epsilon \beta^4). \quad (C.13) $$

At zeroth order the system reads

$$ \partial_{\xi} \eta_0 = \partial_{\xi} u_0, \quad \text{hence } \eta_0 = u_0 \quad \text{together with } \partial_{\xi} \sigma \circ d\mathbf{B}_t^{0,x} = 0. $$

At order $\epsilon$, one has

$$ \begin{aligned} & - \partial_{\xi} u_1 dt + d_{\tau} \eta_0 + \eta_0 \partial_{\xi} \eta_0 dt + \Upsilon^{1/2} \sigma \circ d\mathbf{B}_t^{0,x} \partial_{\xi} \eta_0 + \partial_{\xi} \eta_1 dt + \frac{1}{3} \partial_{\xi \xi \xi}^3 \eta_0 dt = 0, \\ & - \partial_{\xi} \eta_1 dt + d_{\tau} \eta_0 + \Upsilon^{1/2} \sigma \circ d\mathbf{B}_t^{0,x} \partial_{\xi} \eta_0 + 2 \eta_0 \partial_{\xi} \eta_0 dt + \partial_{\xi} u_1 dt + \Upsilon^{1/2} \partial_{\xi} \sigma \circ d\mathbf{B}_t^{1,x} = 0. \end{aligned} $$

Assuming that the noise terms $\sigma \circ d\mathbf{B}_t^{1,x}$ does not depend on space, yet making no assumption on $\sigma \circ d\mathbf{B}_t^{0,x}$, one obtains a stochastic KdV equation with transport noise,

$$ d_{\tau} \eta_0 + \frac{3}{2} \eta_0 \partial_{\xi} \eta_0 dt + \Upsilon^{1/2} \sigma \circ d\mathbf{B}_t^{0,x} \partial_{\xi} \eta_0 + \frac{1}{6} \partial_{\xi \xi \xi}^3 \eta_0 dt = 0. \quad (C.14) $$

Relaxing the spatially constant noise assumption on $\sigma \circ d\mathbf{B}_t^{1,x}$ would add an additive stochastic forcing $\frac{1}{2}\Upsilon^{1/2} \partial_{\xi} \sigma \circ d\mathbf{B}_t^{1,x}$. In such a case, we would get a stochastic KdV equation forced by an additive white noise of the form studied in [41], and for which existence and unicity of solution have been shown in the Sobolev space $H^1(\mathbb{R})$.

Without this additive forcing term, we face a simpler system that boils down to the deterministic one. As a
matter of fact, proceeding to the change of variable suggested by Wadati [42], one has

$$X = \xi - \frac{2}{3}\Upsilon^{1/2} \int_{0}^{t} \boldsymbol{\sigma} \circ d\boldsymbol{B}_{s}^{0,x}, \qquad (C.15)$$

with $\boldsymbol{\sigma}$ being constant over time and space. Now reparametrising $\eta(\xi, \tau)$ as $\eta'(X, t)$, we obtain the unperturbed KdV equation

$$d_t \eta_0' + \frac{3}{2} \eta_0' \partial_x \eta_0' dt + \frac{1}{6} \partial_{xXX}^3 \eta_0' dt = 0, \qquad (C.16)$$

using the chain rules $\partial_{\xi}\eta = \partial_{X}\eta'\partial_{\xi}X = \partial_{X}\eta'$ and $d_{\tau}\eta = \partial_{X}\eta' d_{t}X + d_{t}\eta'$. Such an equation admits a solitary travelling wave solution given by

$$\eta_0'(X, t) = A \operatorname{sech}^2\left(\sqrt{\frac{3}{4}} A \left(X - \left(1 + \frac{1}{2}A\right)t\right)\right). \qquad (C.17)$$

Note that considering a solution with a random phase would lead to the presence of a new additive noise term $\partial_{\varphi}\eta_0 \epsilon f'(B_t)dB_t$. Plus, the noise $\boldsymbol{\sigma} \circ d\boldsymbol{B}_t^{0,x}$ may not be constant in space anymore.

## References

[1] L. Li, B. Deremble, N. Lahaye, 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) (2023) e2022MS003297.

[2] F. L. Tucciarone, E. Mémin, L. Li, Primitive equations under location uncertainty: Analytical description and model development (2023).

[3] E. Mémin, L. Li, N. Lahaye, G. Tissot, B. Chapron, Linear wave solutions of a stochastic shallow water model check for updates, in: Stochastic Transport in Upper Ocean Dynamics II: STUOD 2022 Workshop, London, UK, September 26–29, Vol. 11, Springer Nature, 2023, p. 223.

[4] E. Dinvay, E. Mémin, Hamiltonian formulation of the stochastic surface wave problem, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 478 (2265) (2022) 20220050. doi: 10.1098/rspa.2022.0050.
URL https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2022.0050

[5] E. Mémin, Fluid flow dynamics under location uncertainty, Geophys. & Astro. Fluid Dyn. 108 (2) (2014) 119–146. doi:10.1080/03091929.2013.836190.

[6] J. M. Restrepo, Wave breaking dissipation in the wave-driven ocean circulation, Journal of Physical Oceanography 37 (7) (2007) 1749 – 1763. doi:10.1175/JPO3099.1.
URL https://journals.ametsoc.org/view/journals/phoc/37/7/jpo3099.1.xml

[7] J. M. Restrepo, J. M. Ramírez, J. C. McWilliams, M. Banner, Multiscale momentum flux and diffusion due to whitecapping in wave-current interactions, Journal of Physical Oceanography 41 (5) (2011) 837 – 856. doi:10.1175/2010JP04298.1.
URL https://journals.ametsoc.org/view/journals/phoc/41/5/2010jpo4298.1.xml

[8] P. S. Berloff, J. C. McWilliams, Material transport in oceanic gyres. part ii: Hierarchy of stochastic models, Journal of Physical Oceanography 32 (3) (2002) 797 – 830. doi:10.1175/1520-0485(2002)032<0797:MTIOGP>2.0.CO;2.
URL https://journals.ametsoc.org/view/journals/phoc/32/3/1520-0485_2002_032_0797_mtiogp_2.0.co_2.xml
[9] R. Kraichnan, The structure of isotropic turbulence at very high reynolds numbers, J. of Fluids Mech. 5 (1959) 477-543.

[10] F. Serre, Contribution à l'étude des écoulements permanents et variables dans les canaux, La Houille Blanche (6) (1953) 830-872.

[11] A. E. Green, P. M. Naghdi, A derivation of equations for wave propagation in water of variable depth, Journal of Fluid Mechanics 78 (2) (1976) 237-246.

[12] E. Barthélemy, *Nonlinear shallow water theories for coastal waves*, Surveys in Geophysics 25 (3) (2004) 315-337.
doi:10.1007/s10712-003-1281-7
URL https://doi.org/10.1007/s10712-003-1281-7

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

[14] L. Li, E. Mémin, G. Tissot, Stochastic parameterization with dynamic mode decomposition, in: STUOD Proceedings, Springer Verlag, 2022.

[15] P. Schmid, Dynamic mode decomposition of numerical and experimental data, J. Fluid Mech. 656 (2010) 5-28.

[16] B. Dufée, E. Mémin, D. Crisan, *Stochastic parametrization: An alternative to inflation in ensemble kalman filters*, Quarterly Journal of the Royal Meteorological Society 148 (744) (2022) 1075-1091. doi:https://doi.org/10.1002/qj.4247
URL https://doi.org/10.1002/qj.4247

[17] W. Bauer, P. Chandramouli, B. Chapron, L. Li, E. Mémin, *Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: A stochastic approach*, Journal of Physical Oceanography 50 (4) (01 Apr. 2020) 983-1003. doi:10.1175/JP0-D-19-0164.1
URL https://journals.ametsoc.org/view/journals/phoc/50/4/jpo-d-19-0164.1.xml

[18] H. Kunita, Stochastic flows and stochastic differential equations, Cambridge University Press, 1990.

[19] V. Resseguier, E. Mémin, B. Chapron, Geophysical flows under location uncertainty, Part I Random transport and general models, Geophys. & Astro. Fluid Dyn. 111 (3) (2017) 149-176.

[20] S. Kadri Harouna, E. Mémin, *Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling*, Computers & Fluids 156 (2017) 456-469. doi:10.1016/j.compfluid.2017.08.017.
URL https://hal.inria.fr/hal-01394780

[21] J. Smagorinsky, General circulation experiments with the primitive equation: I. the basic experiment, Monthly Weather Review 91 (1963) 99-165.

[22] M. H. Redi, Oceanic isopycnal mixing by coordinate rotation, Journal of Physical Oceanography 12 (10) (1982) 1154-1158.

[23] A. Craik, S. Leibovich, Rational model for langmuir circulations, J. Fluid Mech. 73 (1976) 401-426.

[24] S. Leibovich, On wave-current interaction theories of langmuir circulations, J. Fluid Mech. 99 (4) (1980) 715-724.

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

[26] E. Barthélemy, *Nonlinear shallow water theories for coastal waves*, Surveys in Geophysics 25 (3) (2004) 315-337. doi:10.1007/s10712-003-1281-7.
URL https://doi.org/10.1007/s10712-003-1281-7
[27] R. Kozlovsky, J. Grobman, The blue garden: coastal infrastructure as ecologically enhanced wave-scapes, Landscape Research (12 2016). doi:10.1080/01426397.2016.1260702.

[28] A. Green, P. Naghdi, A derivation of equations for wave propagation in water of variable depth, Journal of Fluid Mechanics 78 (1976) 237–246.

[29] F. Serre, Contribution à l'étude des écoulements permanents et variables dans les canaux, Houille Blanche 8 (1953) 374–388.

[30] R. Brecht, L. Li, W. Bauer, E. Mémin, Rotating shallow water flow under location uncertainty with a structure-preserving discretization, Journal of advances in modelling earth systems 13 (12) (2021).

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

[32] D. G. Dritschel, M. R. Jalali, On the regularity of the green–naghdi equations for a rotating shallow fluid layer, Journal of Fluid Mechanics 865 (2019) 100–136. doi:10.1017/jfm.2019.47.

[33] G. Vallis, Atmospheric and Oceanic Fluid Dynamics, Cambridge University Press, 2017.

[34] J. Miles, R. Salmon, Weakly dispersive nonlinear gravity waves, Journal of Fluid Mechanics 157 (1985) 519–531.

[35] V. Duchêne, P. Navarro, Waterwaves1d.jl (Jun 2022). URL https://perso.univ-rennes1.fr/vincent.duchene/post/waterwaves1d/

[36] A. Xiao, X. Tang, High strong order stochastic Runge-Kutta methods for Stratonovich stochastic differential equations with scalar noise, Numerical Algorithms 72 (2016) 259–296.

[37] F. Gugole, C. L. E. Franzke, Numerical development and evaluation of an energy conserving conceptual stochastic climate model, Mathematics of Climate and Weather Forecasting 5 (1) (2019) 45–64 [cited 2024-03-21]. doi:doi:10.1515/mcwf-2019-0004.
URL https://doi.org/10.1515/mcwf-2019-0004

[38] G. Pavliotis, A. Stuart, Multiscale methods: averaging and homogenization, Springer Science & Business Media, 2008.

[39] C. Fiorini, P.-M. Boulvard, L. Li, E. Mémin, A two-step numerical scheme in time for surface quasi geostrophic equations under location uncertainty (2023).

[40] R. S. Johnson, A Modern Introduction to the Mathematical Theory of Water Waves, Cambridge Univ. Press, 1997.

[41] A. de Bouard, A. Debussche, On the Stochastic Korteweg de Vries Equation, Journal of Functional Analysis 154 (1998) 215–251.

[42] M. Wadati, Stochastic korteweg-de vries equation, Journal of the Physical Society of Japan 52 (8) (1983) 2642–2648. arXiv:https://doi.org/10.1143/JPSJ.52.2642, doi:10.1143/JPSJ.52.2642.
URL https://doi.org/10.1143/JPSJ.52.2642