# Stochastic Compressible Navier-Stokes Equations Under Location Uncertainty

Gilles Tissot, Étienne Mémin, Quentin Jamet

► To cite this version:

Gilles Tissot, Étienne Mémin, Quentin Jamet. Stochastic Compressible Navier-Stokes Equations Under Location Uncertainty. Stochastic Transport in Upper Ocean Dynamics II, 11, Springer Nature Switzerland, pp.293-319, 2024, Mathematics of Planet Earth, 10.1007/978-3-031-40094-0_14. hal-04356221

HAL Id: hal-04356221

https://inria.hal.science/hal-04356221

Submitted on 10 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Stochastic Compressible Navier-Stokes Equations Under Location Uncertainty

Gilles Tissot, Étienne Mémin, Quentin Jamet

► To cite this version:

Gilles Tissot, Étienne Mémin, Quentin Jamet. Stochastic Compressible Navier-Stokes Equations Under Location Uncertainty. Stochastic Transport in Upper Ocean Dynamics II, 11, Springer Nature Switzerland, pp.293-319, 2024, Mathematics of Planet Earth, 10.1007/978-3-031-40094-0_14 . hal-04356221

HAL Id: hal-04356221

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

Submitted on 10 Apr 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Gilles Tissot, Étienne Mémin, and Quentin Jamet

**Abstract** The aim of this paper is to provide a stochastic version under location uncertainty of the compressible Navier-Stokes equations. To that end, some clarifications of the stochastic Reynolds transport theorem are given when stochastic source terms are present in the right-hand side. We apply this conservation theorem to density, momentum and total energy in order to obtain a transport equation of the primitive variables, i.e. density, velocity and temperature. We show that performing low Mach and Boussinesq approximations to this more general set of equations allows us to recover the known incompressible stochastic Navier-Stokes equations and the stochastic Boussinesq equations, respectively. Finally, we provide some research directions of using this general set of equations in the perspective of relaxing the Boussinesq and hydrostatic assumptions for ocean modelling.

## 1 Introduction

Stochastic modelling under location uncertainty (LU) relies on the decomposition of the displacement of fluid particles into a time-differentiable velocity field, and a highly fluctuating component represented by a Brownian motion. It was proposed in [15] to apply this principle to fluid flows, leading to a stochastic version of the Navier-Stokes equations. On this basis, a similar derivation has been performed for various ocean models, such as Boussinesq models [16], quasi-geostrophic (QG) models [4, 14, 17], surface quasi-geostrophic models (SQG) [18] and shallow water equations [5].

<footnote>G. Tissot (✉) · É. Mémin
INRIA Centre de l'Université de Rennes, IRMAR – UMR CNRS 6625, Rennes, France
e-mail: gilles.tissot@inria.fr</footnote>

<footnote>Q. Jamet
INRIA Centre de l'Université de Rennes, IRMAR – UMR CNRS 6625, Rennes, France
LOPS, Ifremer, Plouzané, France</footnote>
<figure><img src="image_2.png" /></figure>

Gilles Tissot, Étienne Mémin, and Quentin Jamet

**Abstract** The aim of this paper is to provide a stochastic version under location uncertainty of the compressible Navier-Stokes equations. To that end, some clarifications of the stochastic Reynolds transport theorem are given when stochastic source terms are present in the right-hand side. We apply this conservation theorem to density, momentum and total energy in order to obtain a transport equation of the primitive variables, i.e. density, velocity and temperature. We show that performing low Mach and Boussinesq approximations to this more general set of equations allows us to recover the known incompressible stochastic Navier-Stokes equations and the stochastic Boussinesq equations, respectively. Finally, we provide some research directions of using this general set of equations in the perspective of relaxing the Boussinesq and hydrostatic assumptions for ocean modelling.

## 1 Introduction

Stochastic modelling under location uncertainty (LU) relies on the decomposition of the displacement of fluid particles into a time-differentiable velocity field, and a highly fluctuating component represented by a Brownian motion. It was proposed in [15] to apply this principle to fluid flows, leading to a stochastic version of the Navier-Stokes equations. On this basis, a similar derivation has been performed for various ocean models, such as Boussinesq models [16], quasi-geostrophic (QG) models [4, 14, 17], surface quasi-geostrophic models (SQG) [18] and shallow water equations [5].

<footnote>G. Tissot (✉) · É. Mémin
INRIA Centre de l'Université de Rennes, IRMAR – UMR CNRS 6625, Rennes, France
e-mail: gilles.tissot@inria.fr</footnote>

<footnote>Q. Jamet
INRIA Centre de l'Université de Rennes, IRMAR – UMR CNRS 6625, Rennes, France
LOPS, Ifremer, Plouzané, France</footnote>
In the ocean, variations of density, temperature and salinity are of great impor-
tance. In the previously cited models the Boussinesq assumption of small com-
pressibility has been assumed from the start. This is a fair approximation, but
it can become limiting, for instance when radiative transfers heat the ocean
surface. Some research efforts have been performed to account for compressibility
in deterministic oceanic flow models [20, 9, 8] in order to obtain energetically
consistent formulations. Another key aspect is that, Boussinesq models cannot
sustain acoustic waves, which is relevant for two major applications: (i) ocean
acoustics and (ii) numerical simulations of non-Boussinesq models, where pseudo-
compressibility strategies [7] are employed to compute the pressure with explicit
schemes without having to solve an expensive 3D Poisson equation [3]. In addition,
a rigorous development of Boussinesq systems requires to perform the Boussinesq
approximations on the compressible equations [22].

The derivation of a compressible stochastic system cannot be derived from the
incompressible stochastic system, since it corresponds to a generalisation step.
We propose in this paper to start from the classical physical conservation laws
to derive a general stochastic compressible Navier-Stokes system. We verify that
the provided set of equations is consistent with the incompressible stochastic
models previously developed. We will moreover theoretically show some potential
developments enabling to perform a relaxation of the Boussinesq assumption. Such
a procedure will allow us to propose stochastic systems of increasing complexity
lying in between Boussinesq hydrostatic system and a fully compressible flow
dynamics.

The paper is organised as follows. In Sect. 2, we briefly recall the LU formalism
and provide a convenient form of the stochastic Reynolds transport theorem when
the budget of conserved quantities is balanced by external source or flux terms of
stochastic nature. In Sect. 3 we develop the stochastic compressible Navier-Stokes
equations. In Sects. 4 and 5, the low Mach number and Boussinesq approxima-
tions are performed respectively. We verify in these two sections the consistency
with stochastic models previously derived from stochastic isochoric models [16].
This stochastic Boussinesq model is generalised by incorporating thermodynamic
effects. In Sect. 6 these approximations are relaxed and we propose a model which
can be integrated explicitly in time, similarly as in [3]. In Sect. 7 some concluding
remarks are given. In appendix, technical calculation rules and important details to
perform energy budgets are provided.

## 2 Stochastic Reynolds Transport Theorem

The transport of conserved quantities subject to a stochastic transport is described
by the stochastic Reynolds transport theorem (SRTT) introduced in [15]. When
stochastic source terms are involved in the budget, additional covariation terms have
to be taken into account [16]. These terms are usually defined in an implicit manner.
In the ocean, variations of density, temperature and salinity are of great impor-
tance. In the previously cited models the Boussinesq assumption of small com-
pressibility has been assumed from the start. This is a fair approximation, but
it can become limiting, for instance when radiative transfers heat the ocean
surface. Some research efforts have been performed to account for compressibility
in deterministic oceanic flow models [20, 9, 8] in order to obtain energetically
consistent formulations. Another key aspect is that, Boussinesq models cannot
sustain acoustic waves, which is relevant for two major applications: (i) ocean
acoustics and (ii) numerical simulations of non-Boussinesq models, where pseudo-
compressibility strategies [7] are employed to compute the pressure with explicit
schemes without having to solve an expensive 3D Poisson equation [3]. In addition,
a rigorous development of Boussinesq systems requires to perform the Boussinesq
approximations on the compressible equations [22].

The derivation of a compressible stochastic system cannot be derived from the
incompressible stochastic system, since it corresponds to a generalisation step.
We propose in this paper to start from the classical physical conservation laws
to derive a general stochastic compressible Navier-Stokes system. We verify that
the provided set of equations is consistent with the incompressible stochastic
models previously developed. We will moreover theoretically show some potential
developments enabling to perform a relaxation of the Boussinesq assumption. Such
a procedure will allow us to propose stochastic systems of increasing complexity
lying in between Boussinesq hydrostatic system and a fully compressible flow
dynamics.

The paper is organised as follows. In Sect. 2, we briefly recall the LU formalism and provide a convenient form of the stochastic Reynolds transport theorem when the budget of conserved quantities is balanced by external source or flux terms of stochastic nature. In Sect. 3 we develop the stochastic compressible Navier-Stokes equations. In Sects. 4 and 5, the low Mach number and Boussinesq approximations are performed respectively. We verify in these two sections the consistency with stochastic models previously derived from stochastic isochoric models [16]. This stochastic Boussinesq model is generalised by incorporating thermodynamic effects. In Sect. 6 these approximations are relaxed and we propose a model which can be integrated explicitly in time, similarly as in [3]. In Sect. 7 some concluding remarks are given. In appendix, technical calculation rules and important details to perform energy budgets are provided.

2 Stochastic Reynolds Transport Theorem

The transport of conserved quantities subject to a stochastic transport is described
by the stochastic Reynolds transport theorem (SRTT) introduced in [15]. When
stochastic source terms are involved in the budget, additional covariation terms have
to be taken into account [16]. These terms are usually defined in an implicit manner.
In the present section, we briefly present the modelling under location uncertainty, and we rewrite the SRTT in a convenient form for further developments.

In the modelling under location uncertainty [15], the displacement $X(x, t)$ of a particle is written in a differential form as

$$dX(x, t) = u(x, t)dt + \sigma_t dB_t, \quad (1)$$

where $u = (u, v, w)^T$ is a time-differentiable velocity component, and $dB_t$ is the increment of a Brownian motion, whose aim is to model unresolved time-correlated velocity contributions. The correlation operator $\sigma_t$ is an integral operator which involves a spatial convolution in the domain $\Omega$ with a user-defined correlation kernel $\check{\sigma}$, such that

$$(\sigma_t dB_t)^i(x) = \int_{\Omega} \check{\sigma}^{ij}(x, x', t) dB_t^j(x') dx'. \quad (2)$$

Associated with $\sigma_t$, we define the (matrix) variance tensor $\mathbf{a}$ (that corresponds to the one point covariance tensor) such that

$$\mathbf{a}_{ij}(x)dt = \mathbb{E}((\sigma_t dB_t)^i(x) (\sigma_t dB_t)^j(x)). \quad (3)$$

Within this framework, the stochastic transport operator of a scalar quantity $q$ is defined by

$$\mathbb{D}_t q \triangleq d_t q + \left( \underbrace{\left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} + \sigma_t (\nabla \cdot \sigma_t) \right)}_{\mathbf{u}^*} \cdot \nabla \right) q \, dt + (\sigma_t dB_t \cdot \nabla) q - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla q) \, dt, \quad (4)$$

where $\mathbf{u}^*$, called drift velocity, is the resolved velocity corrected by the inhomogeneity and divergence of the noise correlation tensor, respectively. Physical relevance of the drift velocity and the stochastic diffusion $\frac{1}{2}\nabla \cdot (\mathbf{a}\nabla q) \, dt$ has been extensively highlighted in previous studies [e.g. 6, 4].

Variation of $q$ integrated over a transported volume [16] can be written

$$d \int_{\Omega(t)} q \, dx = \int_{\Omega(t)} \left( \mathbb{D}_t q + q \nabla \cdot (\mathbf{u}^* \, dt + \sigma_t dB_t) + \nabla \cdot (\sigma_t \mathbf{h}) \, dt \right) \, dx, \quad (5)$$

with $\mathbf{h}$ defined as follows: when the stochastic transport operator is isolated on the left-hand-side (LHS), $\mathbf{h}$ is associated with the martingale part of the remaining right-hand-side (RHS):

$$\mathbb{D}_t q = f dt + \mathbf{h} \cdot dB_t. \quad (6)$$
In the present section, we briefly present the modelling under location uncertainty, and we rewrite the SRTT in a convenient form for further developments.

In the modelling under location uncertainty [15], the displacement $X(x, t)$ of a particle is written in a differential form as

$$dX(x, t) = u(x, t)dt + \sigma_t dB_t, \quad (1)$$

where $u = (u, v, w)^T$ is a time-differentiable velocity component, and $dB_t$ is the increment of a Brownian motion, whose aim is to model unresolved time-correlated velocity contributions. The correlation operator $\sigma_t$ is an integral operator which involves a spatial convolution in the domain $\Omega$ with a user-defined correlation kernel $\check{\sigma}$, such that

$$(\sigma_t dB_t)^i(x) = \int_{\Omega} \check{\sigma}^{ij}(x, x', t) dB_t^j(x') dx'. \quad (2)$$

Associated with $\sigma_t$, we define the (matrix) variance tensor $\mathbf{a}$ (that corresponds to the one point covariance tensor) such that

$$\mathbf{a}_{ij}(x)dt = \mathbb{E}((\sigma_t dB_t)^i(x) (\sigma_t dB_t)^j(x)). \quad (3)$$

Within this framework, the stochastic transport operator of a scalar quantity $q$ is defined by

$$\mathbb{D}_t q \triangleq d_t q + \left( \underbrace{\left( u - \frac{1}{2} \nabla \cdot \mathbf{a} + \sigma_t (\nabla \cdot \sigma_t) \right)}_{u^*} \cdot \nabla \right) q \, dt + (\sigma_t dB_t \cdot \nabla) q - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla q) \, dt, \quad (4)$$

where $u^*$, called *drift velocity*, is the resolved velocity corrected by the inhomogeneity and divergence of the noise correlation tensor, respectively. Physical relevance of the drift velocity and the stochastic diffusion $\frac{1}{2}\nabla \cdot (\mathbf{a}\nabla q) \, dt$ has been extensively highlighted in previous studies [e.g. 6, 4].

Variation of $q$ integrated over a transported volume [16] can be written

$$d \int_{\Omega(t)} q \, dx = \int_{\Omega(t)} \left( \mathbb{D}_t q + q \nabla \cdot (u^* dt + \sigma_t dB_t) + \nabla \cdot (\sigma_t h) dt \right) dx, \quad (5)$$

with $h$ defined as follows: when the stochastic transport operator is isolated on the left-hand-side (LHS), $h$ is associated with the martingale part of the remaining right-hand-side (RHS):

$$\mathbb{D}_t q = f dt + h \cdot dB_t. \quad (6)$$
Starting from (5), we can assume that some source terms $Q_t dt + Q_\sigma \cdot dB_t$, with a time-differentiable and a martingale contribution respectively, are balancing the budget of $q$ in the control volume such as

$$d \int_{\Omega(t)} q \, dx = \int_{\Omega(t)} (Q_t dt + Q_\sigma \cdot dB_t) \, dx. \qquad (7)$$

These RHS terms correspond to forces (resp. work) when this general expression is associated with the momentum (resp. energy) equation. Dropping the volume integral, we can now identify $h$

$$h \cdot dB_t = -q \nabla \cdot (\sigma_t dB_t) + Q_\sigma \cdot dB_t. \qquad (8)$$

This leads to the explicit expression of the stochastic Reynolds transport theorem

$$d_t q + \nabla \cdot \left( \left( \left( u - \frac{1}{2} \nabla \cdot a \right) dt + \sigma_t dB_t \right) q \right) + \nabla \cdot (\sigma_t Q_\sigma) dt - \frac{1}{2} \nabla \cdot (a \nabla q) dt = Q_t dt + Q_\sigma \cdot dB_t. \qquad (9)$$

The absence of the term $\sigma_t (\nabla \cdot \sigma_t)$ in the modified drift is an important feature of this expression. It has been cancelled (not neglected) by accounting for the term $-q \nabla \cdot (\sigma_t dB_t)$ in Eq. (8). As it will be detailed further, this term will reappear when we will transform the conservative form of the equations to their associated non-conservative form, i.e. writing a transport equation for the primitive variables. For consistency checking, it has been assessed in appendix 7, that the same expression is obtained using a Stratonovich stochastic integral convention.

# 3 Stochastic Compressible Navier-Stokes Equations

To obtain the stochastic compressible Navier-Stokes equations we apply the SRTT equation (9) to the mass, momentum and total energy. This requires at first to properly define the physical variables.

## 3.1 Non-dimensioning

We consider the time $t$, $x = (x, y, z)^T$ the space coordinates of $\Omega$, and $(e_x, e_y, e_z)$ the associated canonical basis. Physical quantities are marked by $\bullet \phi$ and the other quantities are non-dimensional. We adimensionalise by reference conditions (noted $\bullet_{\text{ref}}$), and introduce a reference distance $L_{\text{ref}}$, velocity $u_{\text{ref}}$, density $\rho_{\text{ref}}$, sound speed $c_{\text{ref}}$ as well as viscosity $\mu_{\text{ref}}$. We get
Starting from (5), we can assume that some source terms $Q_t dt + Q_\sigma \cdot dB_t$, with a time-differentiable and a martingale contribution respectively, are balancing the budget of $q$ in the control volume such as

$$d \int_{\Omega(t)} q \, dx = \int_{\Omega(t)} (Q_t dt + Q_\sigma \cdot dB_t) \, dx. \qquad (7)$$

These RHS terms correspond to forces (resp. work) when this general expression is associated with the momentum (resp. energy) equation. Dropping the volume integral, we can now identify $h$

$$h \cdot dB_t = -q \nabla \cdot (\sigma_t dB_t) + Q_\sigma \cdot dB_t. \qquad (8)$$

This leads to the explicit expression of the stochastic Reynolds transport theorem

$$d_t q + \nabla \cdot \left( \left( (u - \frac{1}{2} \nabla \cdot \mathbf{a}) dt + \sigma_t dB_t \right) q \right) + \nabla \cdot (\sigma_t Q_\sigma) dt - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla q) dt = Q_t dt + Q_\sigma \cdot dB_t. \qquad (9)$$

The absence of the term $\sigma_t (\nabla \cdot \sigma_t)$ in the modified drift is an important feature of this expression. It has been cancelled (not neglected) by accounting for the term $-q \nabla \cdot (\sigma_t dB_t)$ in Eq. (8). As it will be detailed further, this term will reappear when we will transform the conservative form of the equations to their associated non-conservative form, i.e. writing a transport equation for the primitive variables. For consistency checking, it has been assessed in appendix 7, that the same expression is obtained using a Stratonovich stochastic integral convention.

# 3 Stochastic Compressible Navier-Stokes Equations

To obtain the stochastic compressible Navier-Stokes equations we apply the SRTT equation (9) to the mass, momentum and total energy. This requires at first to properly define the physical variables.

## 3.1 Non-dimensioning

We consider the time $t$, $x = (x, y, z)^T$ the space coordinates of $\Omega$, and $(e_x, e_y, e_z)$ the associated canonical basis. Physical quantities are marked by $\bullet \phi$ and the other quantities are non-dimensional. We adimensionalise by reference conditions (noted $\bullet_{\text{ref}}$), and introduce a reference distance $L_{\text{ref}}$, velocity $u_{\text{ref}}$, density $\rho_{\text{ref}}$, sound speed $c_{\text{ref}}$ as well as viscosity $\mu_{\text{ref}}$. We get
$$
\begin{align}
\mathbf{x} &= \frac{\mathbf{x}^{\phi}}{L_{\text{ref}}} ; \quad t = \frac{t^{\phi} u_{\text{ref}}}{L_{\text{ref}}} ; \quad \mathbf{u} = \frac{\mathbf{u}^{\phi}}{u_{\text{ref}}} ; \quad c = \frac{c^{\phi}}{u_{\text{ref}}} ; \quad M = \frac{u_{\text{ref}}}{c_{\text{ref}}} ; \nonumber \\
\rho &= \frac{\rho^{\phi}}{\rho_{\text{ref}}} ; \quad \mu = \frac{\mu^{\phi}}{\mu_{\text{ref}}} ; \quad p = \frac{p^{\phi}}{\rho_{\text{ref}} u_{\text{ref}}^{2}} ; \quad T = \frac{T^{\phi} c_{p}^{\phi}}{u_{\text{ref}}^{2}} ; \tag{10} \\
\gamma &= \frac{c_{p}^{\phi}}{c_{v}^{\phi}} ; \quad e = \frac{e^{\phi}}{u_{\text{ref}}^{2}} = \frac{T}{\gamma} ; \quad \mathbf{g} = \frac{\mathbf{g}^{\phi} L_{\text{ref}}}{\rho_{\text{ref}} u_{\text{ref}}^{2}}, \nonumber
\end{align}
$$

with **u** the velocity vector, **c** the speed of sound, **M** the Mach number (i.e. the ratio of typical particle speed to typical sound speed), **ρ** the density, **p** the pressure, **μ** de dynamic viscosity, **T** the temperature, **γ** the heat capacity ratio, (c<sub>p</sub>, c<sub>v</sub>) the heat capacities at constant pressure/volume, **e** the internal energy and **g** = −**g**e<sub>z</sub> the acceleration vector due to gravity. We introduce as well the Reynolds and Prandtl numbers

$$
Re = \frac{\rho_{\text{ref}} u_{\text{ref}} L_{\text{ref}}}{\mu_{\text{ref}}} ; \quad Pr = \frac{c_p^\phi \mu^\phi}{k_T^\phi}, \qquad (11)
$$

with $k_T^\phi$ the thermal conductivity.

3.2 Continuity

Mass conservation ensues upon applying the SRTT on density, i.e. $q = \rho$ with no mass source of any kind:

$$
d_t \rho + \nabla \cdot \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} \right) d_t + \sigma_t d\mathbf{B}_t \right) \rho = \frac{1}{2} \nabla \cdot \left( \mathbf{a} \nabla \rho \right) d_t. \quad (12)
$$

3.3 Momentum

Applying now the SRTT to the momentum ρuᵢ balanced by forces, with uᵢ ∈
{u, v, w}.

$$
\begin{align}
& d_t(\rho u_i) + \nabla \cdot \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t dB_t \right) \rho u_i \nonumber \\
& \quad + \nabla \cdot \left( \sigma_t F_{\sigma}^{\rho u_i} \right) dt \nonumber \\
&= - \frac{\partial p}{\partial x_i} dt - \frac{\partial}{\partial x_i} \left( \rho g \delta_{i,z} \right) \nonumber \\
& \quad + \frac{1}{Re} \frac{\partial \tau_{ij}(\mathbf{u})}{\partial x_j} dt + \frac{1}{Re} \frac{\partial \tau_{ij}(\sigma_t dB_t)}{\partial x_j} dt \nonumber \\
& \quad + \frac{1}{2} \nabla \cdot \left( \mathbf{a} \nabla (\rho u_i) \right) dt, \tag{13}
\end{align}
$$
$$
\begin{align}
\mathbf{x} &= \frac{\mathbf{x}^{\phi}}{L_{\text{ref}}} ; \quad t = \frac{t^{\phi} u_{\text{ref}}}{L_{\text{ref}}} ; \quad \mathbf{u} = \frac{\mathbf{u}^{\phi}}{u_{\text{ref}}} ; \quad c = \frac{c^{\phi}}{u_{\text{ref}}} ; \quad M = \frac{u_{\text{ref}}}{c_{\text{ref}}} ; \nonumber \\
\rho &= \frac{\rho^{\phi}}{\rho_{\text{ref}}} ; \quad \mu = \frac{\mu^{\phi}}{\mu_{\text{ref}}} ; \quad p = \frac{p^{\phi}}{\rho_{\text{ref}} u_{\text{ref}}^{2}} ; \quad T = \frac{T^{\phi} c_{p}^{\phi}}{u_{\text{ref}}^{2}} ; \tag{10} \\
\gamma &= \frac{c_{p}^{\phi}}{c_{v}^{\phi}} ; \quad e = \frac{e^{\phi}}{u_{\text{ref}}^{2}} = \frac{T}{\gamma} ; \quad \mathbf{g} = \frac{\mathbf{g}^{\phi} L_{\text{ref}}}{\rho_{\text{ref}} u_{\text{ref}}^{2}}, \nonumber
\end{align}
$$

with **u** the velocity vector, **c** the speed of sound, **M** the Mach number (i.e. the ratio of typical particle speed to typical sound speed), **ρ** the density, **p** the pressure, **μ** de dynamic viscosity, **T** the temperature, **γ** the heat capacity ratio, (c<sub>p</sub>, c<sub>v</sub>) the heat capacities at constant pressure/volume, **e** the internal energy and **g** = −**g**e<sub>z</sub> the acceleration vector due to gravity. We introduce as well the Reynolds and Prandtl numbers

$$
Re = \frac{\rho_{\text{ref}} u_{\text{ref}} L_{\text{ref}}}{\mu_{\text{ref}}} ; \quad Pr = \frac{c_p^\phi \mu^\phi}{k_T^\phi}, \qquad (11)
$$

with $k_T^\phi$ the thermal conductivity.

3.2 Continuity

Mass conservation ensues upon applying the SRTT on density, i.e. $q = \rho$ with no mass source of any kind:

$$
d_t \rho + \nabla \cdot \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t d\mathbf{B}_t \right) \rho = \frac{1}{2} \nabla \cdot \left( \mathbf{a} \nabla \rho \right) dt. \quad (12)
$$

3.3 Momentum

Applying now the SRTT to the momentum ρuᵢ balanced by forces, with uᵢ ∈
{u, v, w}.

$$
\begin{align}
& d_t (\rho u_i) + \nabla \cdot \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t d\mathbf{B}_t \right) \rho u_i \nonumber \\
& \qquad + \nabla \cdot \left( \sigma_t F_{\sigma}^{\rho u_i} \right) dt \nonumber \\
&= - \frac{\partial p}{\partial x_i} dt - \frac{\partial}{\partial x_i} \left( \rho g \delta_{i,z} \right) \nonumber \\
& \qquad + \frac{1}{Re} \frac{\partial \tau_{ij}(\mathbf{u})}{\partial x_j} dt + \frac{1}{Re} \frac{\partial \tau_{ij}(\sigma_t d\mathbf{B}_t)}{\partial x_j} dt \nonumber \\
& \qquad + \frac{1}{2} \nabla \cdot \left( \mathbf{a} \nabla (\rho u_i) \right) dt, \tag{13}
\end{align}
$$
with $F_{\sigma}^{\rho u_i} \cdot dB_t = -\frac{\partial d p_t^{\sigma}}{\partial x_i} + \frac{1}{Re} \nabla \cdot (\tau_i(\sigma_t dB_t))$. The forces involved here are caused by pressure gradient, viscous stresses $\tau$ and gravity. The pressure gradient is decomposed in a time-differentiable part $p dt$ and a random component $dp_t^{\sigma}$. For sake of generality, we consider the molecular viscosity stress tensor:

$$ \tau(u) = \mu (\nabla u + (\nabla u)^T) + \left(\mu_b - \frac{2}{3}\mu\right) \nabla \cdot u \mathbf{I}, \quad (14) $$

with $\mu_b$ the bulk viscosity. Similarly to the pressure, there is a finite variation friction contribution due to $u dt$ and a martingale contribution due to $\sigma_t dB_t$.

After some manipulations and using the stochastic distributivity rule (70) given in appendix 7, we obtain

$$ \begin{aligned} & \rho \mathbb{D}_t u_i \\ & + \sum_k d_t \left\langle \int_0^t \rho(\sigma_s dB_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} \left( -\frac{\partial d p_s^{\sigma}}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\sigma_s dB_s)}{\partial x_j} \right) \right) \right\rangle \\ & = -\frac{\partial p}{\partial x_i} dt - \frac{\partial d p_t^{\sigma}}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(u)}{\partial x_j} dt + \frac{1}{Re} \frac{\partial \tau_{ij}(\sigma_t dB_t)}{\partial x_j} - \rho g \delta_{i,e_z}. \end{aligned} \quad (15) $$

This expression is very similar to the momentum equation of the incompressible Navier-Stokes equations [15, eq. 41 with incompressibility assumption]. We only have as an additional term the covariation between (the martingale part of) forces and the small scale component $\rho\sigma_t dB_t$. This term, usually difficult to evaluate analytically is generally neglected through a slight variation of the expression of Newton's law in the LU framework, as for instance in [16, Appendix E].

## 3.4 Energy

As in the deterministic framework [23, 2, 13], we now consider conservation of the total energy and deduce a transport equation for the temperature.

### General Formulation

Work of forces and heat fluxes acting on a transported control volume induce variations of total energy E such that:

$$ \begin{aligned} & d_t(\rho E) + \nabla \cdot \left( \left( (u - \frac{1}{2} \nabla \cdot \mathbf{a}) dt + \sigma_t dB_t \right) \rho E \right) + \nabla \cdot (\sigma_t F_{\sigma}^{\rho E}) \\ & = \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla(\rho E)) dt + dW - \nabla \cdot (dq), \end{aligned} \quad (16) $$
with $F_{\sigma}^{\rho u_i} \cdot dB_t = -\frac{\partial d p_t^{\sigma}}{\partial x_i} + \frac{1}{Re} \nabla \cdot (\tau_i (\sigma_t dB_t))$. The forces involved here are caused by pressure gradient, viscous stresses $\tau$ and gravity. The pressure gradient is decomposed in a time-differentiable part $p dt$ and a random component $dp_t^{\sigma}$. For sake of generality, we consider the molecular viscosity stress tensor:

$$ \tau(u) = \mu (\nabla u + (\nabla u)^T) + \left(\mu_b - \frac{2}{3}\mu\right) \nabla \cdot u \mathbf{I}, \quad (14) $$

with $\mu_b$ the bulk viscosity. Similarly to the pressure, there is a finite variation friction contribution due to $u dt$ and a martingale contribution due to $\sigma_t dB_t$.

After some manipulations and using the stochastic distributivity rule (70) given in appendix 7, we obtain

$$ \begin{aligned} & \rho \mathbb{D}_t u_i \\ & \quad + \sum_k d_t \left\langle \int_0^t \rho (\sigma_s dB_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} \left( -\frac{\partial d p_s^{\sigma}}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\sigma_s dB_s)}{\partial x_j} \right) \right) \right\rangle \\ & = -\frac{\partial p}{\partial x_i} dt - \frac{\partial d p_t^{\sigma}}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(u)}{\partial x_j} dt + \frac{1}{Re} \frac{\partial \tau_{ij}(\sigma_t dB_t)}{\partial x_j} - \rho g \delta_{i,e_z}. \end{aligned} \quad (15) $$

This expression is very similar to the momentum equation of the incompressible Navier-Stokes equations [15, eq. 41 with incompressibility assumption]. We only have as an additional term the covariation between (the martingale part of) forces and the small scale component $\rho \sigma_t dB_t$. This term, usually difficult to evaluate analytically is generally neglected through a slight variation of the expression of Newton's law in the LU framework, as for instance in [16, Appendix E].

## 3.4 Energy

As in the deterministic framework [23, 2, 13], we now consider conservation of the total energy and deduce a transport equation for the temperature.

### General Formulation

Work of forces and heat fluxes acting on a transported control volume induce variations of total energy E such that:

$$ \begin{aligned} & d_t(\rho E) + \nabla \cdot \left( \left( (u - \frac{1}{2} \nabla \cdot \mathbf{a}) dt + \sigma_t dB_t \right) \rho E \right) + \nabla \cdot (\sigma_t F_{\sigma}^{\rho E}) \\ & = \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla (\rho E)) dt + dW - \nabla \cdot (dq), \end{aligned} \quad (16) $$
with $dW$ and $dq$ the elementary work of the forces and heat fluxes detailed later. The martingale part of these RHS terms is written $F_{\sigma}^{\rho E} \cdot dB_t$.

Using (70) and the continuity equation (12), we obtain

$$ \rho \mathbb{D}_t(E) + \sum_k d_t \left\langle \int_0^t \rho(\sigma_s dB_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} F_{\sigma}^{\rho E} \cdot dB_s \right) \right\rangle = dW - \nabla \cdot (dq). \quad (17) $$

## Definition of the Energy

At this point, the form of the total energy has to be specified. It is strongly related to the physical mechanisms at play. In the present study, we consider the total energy $\rho E = \rho(e + \frac{1}{2}\|u\|^2 + gz)$, as the sum of internal energy $e = \frac{T}{\gamma}$, kinetic energy and potential energy due to gravity. We do not consider the energy of the Brownian motion since it is possibly infinite.

## Definition of the Work of Forces and Heat Fluxes

The work of the time-differentiable pressure represents how pressure is working with the displacement of the control surface. The expression can be obtained by integrating the force multiplied by the surface displacement over a transported control volume and applying Green's formulae. The procedure is similar to the deterministic framework, with the additional implication of the drift velocity, as demonstrated in appendix 7. We have for the pressure work:

$$ \begin{aligned} \int_{\Omega(t)} dW_p \, dx &= \int_{\delta\Omega(t)} (-p \, n \, dS) \cdot (u^* \, dt + \sigma_t \, dB_t) \\ &= - \int_{\Omega(t)} \nabla \cdot (p(u^* \, dt + \sigma_t \, dB_t)) \, dx. \end{aligned} \quad (18) $$

The minus sign comes from the outward normal $n$ convention. We can then identify

$$ dW_p = -\nabla \cdot (p(u^* dt + \sigma_t dB_t)). \quad (19) $$

In the same way, the viscous stress of the resolved component can be written

$$ dW_{\tau} = \frac{1}{Re} \nabla \cdot (\tau(u)(u^* dt + \sigma_t dB_t)). \quad (20) $$

Following Appendix 7, we take as well into account the work of the random pressure:

$$ dW_{rp} = -\nabla \cdot (u^* d p_t^{\sigma}), \quad (21) $$

and the work of the random viscous stress
with $dW$ and $dq$ the elementary work of the forces and heat fluxes detailed later. The martingale part of these RHS terms is written $F_{\sigma}^{\rho E} \cdot dB_t$.

Using (70) and the continuity equation (12), we obtain

$$ \rho \mathbb{D}_t(E) + \sum_k d_t \left\langle \int_0^t \rho(\sigma_s dB_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} F_{\sigma}^{\rho E} \cdot dB_s \right) \right\rangle = dW - \nabla \cdot (dq). \quad (17) $$

## Definition of the Energy

At this point, the form of the total energy has to be specified. It is strongly related to the physical mechanisms at play. In the present study, we consider the total energy $\rho E = \rho(e + \frac{1}{2}\|u\|^2 + gz)$, as the sum of internal energy $e = \frac{T}{\gamma}$, kinetic energy and potential energy due to gravity. We do not consider the energy of the Brownian motion since it is possibly infinite.

## Definition of the Work of Forces and Heat Fluxes

The work of the time-differentiable pressure represents how pressure is working with the displacement of the control surface. The expression can be obtained by integrating the force multiplied by the surface displacement over a transported control volume and applying Green's formulae. The procedure is similar to the deterministic framework, with the additional implication of the drift velocity, as demonstrated in appendix 7. We have for the pressure work:

$$ \begin{aligned} \int_{\Omega(t)} dW_p \, dx &= \int_{\delta\Omega(t)} (-p \, n \, dS) \cdot (u^* dt + \sigma_t dB_t) \\ &= - \int_{\Omega(t)} \nabla \cdot (p(u^* dt + \sigma_t dB_t)) \, dx. \end{aligned} \quad (18) $$

The minus sign comes from the outward normal $n$ convention. We can then identify

$$ dW_p = -\nabla \cdot (p(u^* dt + \sigma_t dB_t)). \quad (19) $$

In the same way, the viscous stress of the resolved component can be written

$$ dW_{\tau} = \frac{1}{Re} \nabla \cdot (\tau(u)(u^* dt + \sigma_t dB_t)). \quad (20) $$

Following Appendix 7, we take as well into account the work of the random pressure:

$$ dW_{rp} = -\nabla \cdot (u^* d p_t^{\sigma}), \quad (21) $$

and the work of the random viscous stress
$$dW_{r\tau} = \frac{1}{Re} \nabla \cdot (\tau(\sigma_t d\mathbf{B}_t) \mathbf{u}^*). \qquad (22)$$

As rigorously detailed in appendix 7, we do not consider work of random forces associated with $\sigma_t d\mathbf{B}_t$, since such a work would be highly irregular (in time) and should be in balance with variations of kinetic energy of $\sigma_t d\mathbf{B}_t$, which is possibly infinite and not described by the present model.

There is no work contribution of gravity on total energy, since the gain in kinetic energy directly associated with the gravity force is compensated by the loss in potential energy.

Finally, we obtain the thermal conductivity by expressing the thermal fluxes by the Fourier law $d\mathbf{q} = -\frac{1}{RePr} \nabla T dt$.

## Transport Equation of Temperature

By replacing the energy by the contributions of internal, kinetic and potential energy, and by subtracting the contribution of the kinetic energy using the momentum equation (15) and the distributivity rule (70), we obtain the transport equation for the temperature

$$ \begin{aligned} & \frac{\rho}{\gamma} \mathbb{D}_t T + \underbrace{\sum_k d_t \left( \int_0^t \rho(\sigma_s d\mathbf{B}_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{\gamma}{\rho} \mathbf{F}_{\sigma}^T \cdot \boldsymbol{\sigma}_s d\mathbf{B}_s \right) \right)}_{Q_T} \\ & + \underbrace{\sum_i \frac{\rho}{2} d_t \left( \int_0^t \mathbf{F}_{\sigma}^{u_i} \cdot \boldsymbol{\sigma}_s d\mathbf{B}_s, \int_0^t \mathbf{F}_{\sigma}^{u_i} \cdot \boldsymbol{\sigma}_s d\mathbf{B}_s \right)}_{Q_u} \\ & = \underbrace{-p \nabla \cdot (\mathbf{u}^* dt + \boldsymbol{\sigma}_t d\mathbf{B}_t)}_{P_t} - \underbrace{d p_t^{\sigma} \nabla \cdot \mathbf{u}^*}_{P_{\sigma}} \\ & + \underbrace{\frac{1}{Re} \tau(\mathbf{u}) : \nabla (\mathbf{u}^* dt + \boldsymbol{\sigma}_t d\mathbf{B}_t)}_{V_t} + \underbrace{\frac{1}{Re} \tau(\boldsymbol{\sigma}_t d\mathbf{B}_t) : \nabla \mathbf{u}^*}_{V_{\sigma}} \qquad (23) \\ & + \underbrace{((\mathbf{u}^* - \mathbf{u}) dt + \boldsymbol{\sigma}_t d\mathbf{B}_t) \cdot \left( -\nabla p + \frac{1}{Re} \nabla \cdot \boldsymbol{\tau}(\mathbf{u}) + \rho g \right)}_{D_t} \\ & + \underbrace{(\mathbf{u}^* - \mathbf{u}) \cdot \left( -\nabla d p_t^{\sigma} + \frac{1}{Re} \nabla \cdot \boldsymbol{\tau}(\boldsymbol{\sigma}_t d\mathbf{B}_t) \right)}_{D_{\sigma}} \\ & + \frac{1}{RePr} \nabla \cdot (\nabla T) dt. \end{aligned} $$
$$dW_{r\tau} = \frac{1}{Re} \nabla \cdot (\tau(\sigma_t d\mathbf{B}_t)\mathbf{u}^*). \quad (22)$$

As rigorously detailed in appendix 7, we do not consider work of random forces associated with $\sigma_t d\mathbf{B}_t$, since such a work would be highly irregular (in time) and should be in balance with variations of kinetic energy of $\sigma_t d\mathbf{B}_t$, which is possibly infinite and not described by the present model.

There is no work contribution of gravity on total energy, since the gain in kinetic energy directly associated with the gravity force is compensated by the loss in potential energy.

Finally, we obtain the thermal conductivity by expressing the thermal fluxes by the Fourier law $d\mathbf{q} = -\frac{1}{RePr} \nabla T dt$.

## Transport Equation of Temperature

By replacing the energy by the contributions of internal, kinetic and potential energy, and by subtracting the contribution of the kinetic energy using the momentum equation (15) and the distributivity rule (70), we obtain the transport equation for the temperature

$$ \begin{aligned} & \frac{\rho}{\gamma} \mathbb{D}_t T + \underbrace{\sum_k d_t \left( \int_0^t \rho(\sigma_s d\mathbf{B}_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{\gamma}{\rho} \mathbf{F}_{\sigma}^T \cdot \boldsymbol{\sigma}_s d\mathbf{B}_s \right) \right)}_{Q_T} \\ & + \underbrace{\sum_i \frac{\rho}{2} d_t \left( \int_0^t \mathbf{F}_{\sigma}^{u_i} \cdot \boldsymbol{\sigma}_s d\mathbf{B}_s, \int_0^t \mathbf{F}_{\sigma}^{u_i} \cdot \boldsymbol{\sigma}_s d\mathbf{B}_s \right)}_{Q_u} \\ & = \underbrace{-p \nabla \cdot (\mathbf{u}^* dt + \boldsymbol{\sigma}_t d\mathbf{B}_t)}_{P_t} - \underbrace{d p_t^{\sigma} \nabla \cdot \mathbf{u}^*}_{P_{\sigma}} \\ & + \underbrace{\frac{1}{Re} \tau(\mathbf{u}) : \nabla (\mathbf{u}^* dt + \boldsymbol{\sigma}_t d\mathbf{B}_t)}_{V_t} + \underbrace{\frac{1}{Re} \tau(\boldsymbol{\sigma}_t d\mathbf{B}_t) : \nabla \mathbf{u}^*}_{V_{\sigma}} \quad (23) \\ & + \underbrace{((\mathbf{u}^* - \mathbf{u}) dt + \boldsymbol{\sigma}_t d\mathbf{B}_t) \cdot \left( -\nabla p + \frac{1}{Re} \nabla \cdot \boldsymbol{\tau}(\mathbf{u}) + \rho g \right)}_{D_t} \\ & + \underbrace{(\mathbf{u}^* - \mathbf{u}) \cdot \left( -\nabla d p_t^{\sigma} + \frac{1}{Re} \nabla \cdot \boldsymbol{\tau}(\boldsymbol{\sigma}_t d\mathbf{B}_t) \right)}_{D_{\sigma}} \\ & + \frac{1}{RePr} \nabla \cdot (\nabla T) dt. \end{aligned} $$
with $F_{\sigma}^{u_i} = \frac{1}{\rho} F_{\sigma}^{\rho u_i}$, and $\frac{\gamma}{\rho} F_{\sigma}^{T} \cdot dB_t$ the sum of all martingale terms of the RHS of Eq. (23). In Eq. (23), we recover the terms present in the deterministic framework, but considering the stochastic transport operator instead of the deterministic transport operator. Nevertheless, some covariation terms are now arising. In particular the term $Q_T$ is induced by the random work of the forces, and the term $Q_u$ is induced by the increase of kinetic energy through covariations of the forces in the momentum equation. On the RHS, we remark that the drift velocity is involved in the work of the time-differentiable pressure (in $P_t$) and random pressure (in $P_{\sigma}$), consistently with Appendix 7. The terms $V_t$ and $V_{\sigma}$ are smooth in time and random viscous stresses, respectively. In addition, the terms $D_t$ and $D_{\sigma}$ correspond to works caused by the alignment between the drift and random velocities with the pressure gradient and viscous forces. We call them drift works. Focusing on $-(u^* - u) \cdot \nabla p$, we interpret this drift work to be related to baropycnal work [1], present in the compressible large-eddy simulation framework. Indeed in standard compressible LES, baropycnal work corresponds to a contribution caused by the alignment between the large scale pressure gradient and the Reynolds stresses induced by product between the small scales contributions of $\rho$ and $u$ (i.e. $\frac{1}{\rho}\nabla p \cdot \overline{\rho'u'}$, with ' denoting here small scale components and · large scale filtering). This Reynolds stress $\frac{1}{\rho}\overline{\rho'u'}$ has the dimension of a velocity. In our case, the interpretation of the effective displacement associated with this work is directly the drift velocity over $dt$. Similar interpretations can be made for the other drift work terms, associated with viscous stresses and random variables. The presence of gravity in the drift work $D_t$ shows that in the vertical direction the time-differentiable drift work is of the form $(w^* - w)(\frac{\partial p}{\partial z} - \rho g)$, and we see appearing the vertical small-scale mass flux times the buoyancy $(w^* - w)(\rho_0 b)$, plus non-hydrostatic pressure effects. It can be noticed that for a divergence-free homogeneous noise (for which the variance tensor is constant in space) the drift work is null as $u^* - u$ cancels.

## 3.5 Equation of State

In order to close the system, we have to specify the equation of state. We keep generality and write the equation of state formally as follows

$$p = f(\rho, T). \qquad (24)$$

As in the deterministic framework, since we have an evolution equation of density and temperature, the pressure can be determined explicitly, at the price of a Courant-Friedrichs-Lewi (CFL) condition constrained by the speed of sound.

The random pressure can be identified by differentiating the equation of state. Indeed, we have an explicit evolution equation of the pressure, which can be expressed through Itô formulae (the equation of state $f$ being deterministic—i.e the state map does not depend on the random events) as
with $F_{\sigma}^{u_i} = \frac{1}{\rho} F_{\sigma}^{\rho u_i}$, and $\frac{\gamma}{\rho} F_{\sigma}^{T} \cdot dB_t$ the sum of all martingale terms of the RHS of Eq. (23). In Eq. (23), we recover the terms present in the deterministic framework, but considering the stochastic transport operator instead of the deterministic transport operator. Nevertheless, some covariation terms are now arising. In particular the term $Q_T$ is induced by the random work of the forces, and the term $Q_u$ is induced by the increase of kinetic energy through covariations of the forces in the momentum equation. On the RHS, we remark that the drift velocity is involved in the work of the time-differentiable pressure (in $P_t$) and random pressure (in $P_{\sigma}$), consistently with Appendix 7. The terms $V_t$ and $V_{\sigma}$ are smooth in time and random viscous stresses, respectively. In addition, the terms $D_t$ and $D_{\sigma}$ correspond to works caused by the alignment between the drift and random velocities with the pressure gradient and viscous forces. We call them drift works. Focusing on $-(u^* - u) \cdot \nabla p$, we interpret this drift work to be related to baropycnal work [1], present in the compressible large-eddy simulation framework. Indeed in standard compressible LES, baropycnal work corresponds to a contribution caused by the alignment between the large scale pressure gradient and the Reynolds stresses induced by product between the small scales contributions of $\rho$ and $u$ (i.e. $\frac{1}{\rho}\nabla p \cdot \overline{\rho'u'}$, with ' denoting here small scale components and · large scale filtering). This Reynolds stress $\frac{1}{\rho}\overline{\rho'u'}$ has the dimension of a velocity. In our case, the interpretation of the effective displacement associated with this work is directly the drift velocity over $dt$. Similar interpretations can be made for the other drift work terms, associated with viscous stresses and random variables. The presence of gravity in the drift work $D_t$ shows that in the vertical direction the time-differentiable drift work is of the form $(w^* - w)(\frac{\partial p}{\partial z} - \rho g)$, and we see appearing the vertical small-scale mass flux times the buoyancy $(w^* - w)(\rho_0 b)$, plus non-hydrostatic pressure effects. It can be noticed that for a divergence-free homogeneous noise (for which the variance tensor is constant in space) the drift work is null as $u^* - u$ cancels.

## 3.5 Equation of State

In order to close the system, we have to specify the equation of state. We keep generality and write the equation of state formally as follows

$$p = f(\rho, T). \qquad (24)$$

As in the deterministic framework, since we have an evolution equation of density and temperature, the pressure can be determined explicitly, at the price of a Courant-Friedrichs-Lewi (CFL) condition constrained by the speed of sound.

The random pressure can be identified by differentiating the equation of state. Indeed, we have an explicit evolution equation of the pressure, which can be expressed through Itô formulae (the equation of state $f$ being deterministic—i.e the state map does not depend on the random events) as
$$ \begin{aligned} d_t p &= \frac{\partial f}{\partial \rho} d_t \rho + \frac{\partial f}{\partial T} d_t T + \frac{1}{2} \frac{\partial^2 f}{\partial \rho^2} d_t \langle \rho, \rho \rangle + \frac{1}{2} \frac{\partial^2 f}{\partial T^2} d_t \langle T, T \rangle + \frac{\partial^2 f}{\partial \rho \partial T} d_t \langle \rho, T \rangle \\ &= \frac{\partial \tilde{p}}{\partial t} dt + \frac{d p_t^\sigma}{\tau}, \end{aligned} \quad (25) $$

where $\tilde{p}$ is the time-differentiable part of the pressure which contains, among other things, all covariation terms. The martingale part of $d_t p$ is $d p_t^\sigma / \tau$, with $\tau$ a decorrelation time. This decorrelation time represents the typical time during which the random pressure acts in a coherent manner to produce a change of momentum. It is assumed to be the same decorrelation time than the one classically introduced [e.g. 12, 4] to relate in practice the definition of the variance tensor to velocity fluctuations variance: (i.e. $\mathbf{a} = \tau \mathbb{E}(\mathbf{u}'\mathbf{u}'^T)$). The term $d p^\sigma$ is identified from (25) to be the random pressure acting on the momentum equation.

If we assume that the random pressure ensues from an isentropic process, i.e. of acoustic nature, we can write

$$ d_t p = \left. \frac{\partial p}{\partial \rho} \right|_s d_t \rho = c^2 d_t \rho, \quad (26) $$

with $c$ the speed of sound and $s$ the entropy. We can then identify from (12)

$$ d p_t^\sigma = -\tau c^2 \nabla \cdot (\rho \sigma_t d \mathbf{B}_t). \quad (27) $$

It can be remarked that this expression is consistent with Eq. (25) under the isentropic transformation assumption.

For oceanic flows, the equation of state is often expressed in terms of density rather in pressure. A specific treatment adapted to oceanic flows is detailed in Sect. 6.

## 4 Low Mach Approximation

To perform the low Mach approximation, we follow the same steps as [11], but applied to the compressible stochastic Navier-Stokes equations. With our non-dimensioning, we have at infinity for isentropic transformations,

$$ \left. \frac{\partial p}{\partial \rho} \right|_s = c_{\text{ref}}^2 = \frac{1}{M^2}. \quad (28) $$

This suggests for small $M$ the following asymptotic expansion
$$ 
\begin{aligned}
d_t p &= \frac{\partial f}{\partial \rho} d_t \rho + \frac{\partial f}{\partial T} d_t T + \frac{1}{2} \frac{\partial^2 f}{\partial \rho^2} d_t \langle \rho, \rho \rangle + \frac{1}{2} \frac{\partial^2 f}{\partial T^2} d_t \langle T, T \rangle + \frac{\partial^2 f}{\partial \rho \partial T} d_t \langle \rho, T \rangle \\
&= \frac{\partial \tilde{p}}{\partial t} dt + \frac{d p_t^\sigma}{\tau},
\end{aligned} 
\quad (25) $$

where $\tilde{p}$ is the time-differentiable part of the pressure which contains, among other things, all covariation terms. The martingale part of $d_t p$ is $d p_t^\sigma / \tau$, with $\tau$ a decorrelation time. This decorrelation time represents the typical time during which the random pressure acts in a coherent manner to produce a change of momentum. It is assumed to be the same decorrelation time than the one classically introduced [e.g. 12, 4] to relate in practice the definition of the variance tensor to velocity fluctuations variance: (i.e. $\mathbf{a} = \tau \mathbb{E}(\mathbf{u}'\mathbf{u}'^T)$). The term $d p^\sigma$ is identified from (25) to be the random pressure acting on the momentum equation.

If we assume that the random pressure ensues from an isentropic process, i.e. of acoustic nature, we can write

$$ d_t p = \left. \frac{\partial p}{\partial \rho} \right|_s d_t \rho = c^2 d_t \rho, \quad (26) $$

with $c$ the speed of sound and $s$ the entropy. We can then identify from (12)

$$ d p_t^\sigma = -\tau c^2 \nabla \cdot (\rho \sigma_t d \mathbf{B}_t). \quad (27) $$

It can be remarked that this expression is consistent with Eq. (25) under the isentropic transformation assumption.

For oceanic flows, the equation of state is often expressed in terms of density rather in pressure. A specific treatment adapted to oceanic flows is detailed in Sect. 6.

## 4 Low Mach Approximation

To perform the low Mach approximation, we follow the same steps as [11], but applied to the compressible stochastic Navier-Stokes equations. With our non-dimensioning, we have at infinity for isentropic transformations,

$$ \left. \frac{\partial p}{\partial \rho} \right|_s = c_{\text{ref}}^2 = \frac{1}{M^2}. \quad (28) $$

This suggests for small $M$ the following asymptotic expansion
$$
\begin{align}
& \rho = \rho_0 + M^2 \rho_1 + o(M^2), \nonumber \\
& \mathbf{u} = \mathbf{u}_0 + o(1), \tag{29} \\
& T = \frac{1}{M^2} T_0 + T_1 + o(1). \nonumber
\end{align}
$$

and $p = O(\frac{\rho_1}{M^2}) = O(1)$. Similarly, the random pressure $dp_t^\sigma$ follows the same scaling as the time-differentiable pressure.

Collecting $O(1)$ and $O(M^2)$ terms in the continuity equation, we obtain respectively

$$
\nabla \cdot (\boldsymbol{u}_0^* \mathrm{d}t + \boldsymbol{\sigma}_t \mathrm{d}\boldsymbol{B}_t) = 0 \quad ; \quad \mathbb{D}_t \rho_1 = 0. \tag{30}
$$

In the momentum equation (15), the order of magnitude of the covariation
term can be determined by integrating over the domain, using distributivity of the
divergence and performing an integration by parts:

$$
\begin{align}
& \int_{\Omega} \sum_k d_t \left\langle \int_0^t \rho(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \rho \left(-\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j}\right) \right) \right\rangle \mathrm{d}\mathbf{x} \nonumber \\
&= \int_{\Omega} \sum_k \frac{\partial}{\partial x_k} d_t \left\langle \int_0^t \rho(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)^k, \int_0^t \left( \rho \left(-\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j}\right) \right) \right\rangle \mathrm{d}\mathbf{x} \nonumber \\
& \quad - \int_{\Omega} d_t \left\langle \int_0^t \nabla \cdot (\rho \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s), \int_0^t \left( \rho \left(-\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j}\right) \right) \right\rangle \mathrm{d}\mathbf{x} \nonumber \\
&= \int_{\delta\Omega} d_t \left\langle \int_0^t \rho \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \boldsymbol{n}, \int_0^t \left( \rho \left(-\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j}\right) \right) \right\rangle \mathrm{d}S \nonumber \\
& \quad - \int_{\Omega} d_t \left\langle \int_0^t \nabla \cdot (\rho \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s), \int_0^t \left( \rho \left(-\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j}\right) \right) \right\rangle \mathrm{d}\mathbf{x} \nonumber \\
&= \int_{\Omega} d_t \left\langle \int_0^t \underbrace{d_t \rho}_{\mathcal{O}(M^2)}, \int_0^t \underbrace{\frac{1}{\rho} \mathbf{F}_{\sigma}^{\rho u_i} \cdot \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s}_{\mathcal{O}(1)} \right\rangle \mathrm{d}\mathbf{x} = \mathcal{O}(M^2), \tag{31}
\end{align}
$$

where suitable boundary conditions at $\delta\Omega$ (e.g. Dirichlet boundary conditions (no random inflow velocity) or zero normal stress (outflow boundary conditions)), have been applied to insure the first surface term vanishes.

By neglecting the order $O(M^2)$ terms, we obtain then finally the incompressible
Navier-Stokes presented in [15] under the incompressibility assumption

$$
\begin{align}
& \rho_0 \mathbb{D}_t \mathbf{u} = -\nabla p \, dt - \nabla d p_t^\sigma + \frac{1}{Re} \nabla \cdot (\boldsymbol{\tau}(\mathbf{u})) \, dt + \frac{1}{Re} \nabla \cdot (\boldsymbol{\tau}(\boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t)) + \rho g. \nonumber \\
& \nabla \cdot \mathbf{u}^* = 0 \quad ; \quad \nabla \cdot (\boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t) = 0. \tag{32}
\end{align}
$$
$$ 
\begin{aligned}
& \rho = \rho_0 + M^2 \rho_1 + o(M^2), \\
& \mathbf{u} = \mathbf{u}_0 + o(1), \\
& T = \frac{1}{M^2} T_0 + T_1 + o(1).
\end{aligned} 
\quad (29) $$

and $p = O(\frac{\rho_1}{M^2}) = O(1)$. Similarly, the random pressure $dp_t^\sigma$ follows the same scaling as the time-differentiable pressure.

Collecting $O(1)$ and $O(M^2)$ terms in the continuity equation, we obtain respectively

$$ \nabla \cdot (\mathbf{u}_0^* \mathrm{d}t + \boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t) = 0 \quad ; \quad \mathbb{D}_t \rho_1 = 0. \quad (30) $$

In the momentum equation (15), the order of magnitude of the covariation term can be determined by integrating over the domain, using distributivity of the divergence and performing an integration by parts:

$$ 
\begin{aligned}
& \int_{\Omega} \sum_k d_t \left\langle \int_0^t \rho(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \rho \left( -\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j} \right) \right) \right\rangle \mathrm{d}\mathbf{x} \\
&= \int_{\Omega} \sum_k \frac{\partial}{\partial x_k} d_t \left\langle \int_0^t \rho(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)^k, \int_0^t \left( \rho \left( -\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j} \right) \right) \right\rangle \mathrm{d}\mathbf{x} \\
&\quad - \int_{\Omega} d_t \left\langle \int_0^t \nabla \cdot (\rho \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s), \int_0^t \left( \rho \left( -\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j} \right) \right) \right\rangle \mathrm{d}\mathbf{x} \\
&= \int_{\delta\Omega} d_t \left\langle \int_0^t \rho \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s \cdot \mathbf{n}, \int_0^t \left( \rho \left( -\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j} \right) \right) \right\rangle \mathrm{d}S \\
&\quad - \int_{\Omega} d_t \left\langle \int_0^t \nabla \cdot (\rho \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s), \int_0^t \left( \rho \left( -\frac{\partial \mathrm{d}p_s^\sigma}{\partial x_i} + \frac{1}{Re} \frac{\partial \tau_{ij}(\boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s)}{\partial x_j} \right) \right) \right\rangle \mathrm{d}\mathbf{x} \\
&= \int_{\Omega} d_t \left\langle \int_0^t \underbrace{d_t \rho}_{\mathcal{O}(M^2)}, \int_0^t \underbrace{\frac{1}{\rho} \mathbf{F}_\sigma^{\rho u_i} \cdot \boldsymbol{\sigma}_s \mathrm{d}\mathbf{B}_s}_{\mathcal{O}(1)} \right\rangle \mathrm{d}\mathbf{x} = O(M^2),
\end{aligned} 
\quad (31) $$

where suitable boundary conditions at $\delta\Omega$ (e.g. Dirichlet boundary conditions (no random inflow velocity) or zero normal stress (outflow boundary conditions)), have been applied to insure the first surface term vanishes.

By neglecting the order $O(M^2)$ terms, we obtain then finally the incompressible Navier-Stokes presented in [15] under the incompressibility assumption

$$ 
\begin{aligned}
\rho_0 \mathbb{D}_t \mathbf{u} &= -\nabla p \mathrm{d}t - \nabla \mathrm{d}p_t^\sigma + \frac{1}{Re} \nabla \cdot (\boldsymbol{\tau}(\mathbf{u})) \mathrm{d}t + \frac{1}{Re} \nabla \cdot (\boldsymbol{\tau}(\boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t)) + \rho g. \\
\nabla \cdot \mathbf{u}^* &= 0 \quad ; \quad \nabla \cdot (\boldsymbol{\sigma}_t \mathrm{d}\mathbf{B}_t) = 0.
\end{aligned} 
\quad (32) $$
# 5 Boussinesq-Hydrostatic Approximation

In this section, starting from the stochastic compressible Navier-Stokes equations, we perform the Boussinesq approximation by considering small density fluctuations. These fluctuations are neglected, when they are not multiplied by gravity $g$, which leads to the classical definition of the buoyancy. We perform as well the hydrostatic approximation through the classical aspect ratio scaling $D = H/L_{ref} \ll 1$, with $H$ the water depth. For simplicity, we do not consider a rotating frame. Coriolis correction could be straightforwardly considered as in [21]. The vertical coordinate $z \in [-H, \eta]$ is bounded by the bottom and the free surface.

## Density

The density is decomposed through the following asymptotic expansion

$$ \rho = \rho_0 + \epsilon\rho_1(z) + \epsilon\rho_2(x, y, z, t) + o(\epsilon), \quad (33) $$

with $\rho_1(z)$ the time-averaged stratification term, and $\epsilon \ll 1$ and we do not need to assume that $\rho_1 > \rho_2$. We obtain hence

$$ \nabla \cdot \mathbf{u}^* = 0 \quad ; \quad \nabla \cdot \boldsymbol{\sigma}_t \, d\mathbf{B}_t = 0 \quad ; \quad \mathbb{D}_t (\rho_1 + \rho_2) = 0. \quad (34) $$

The drift velocity and the noise are divergence free. Density perturbations undergo a stochastic transport by the flow. We remark that since $\nabla \cdot \boldsymbol{\sigma}_t d\mathbf{B}_t = 0$, then the transport operator $\mathbb{D}_t(\cdot)$ can be directly used.

The terms of order $\epsilon$ of Eq. (34) can be expressed in terms of buoyancy $b = -\epsilon g\rho_2/\rho_0$:

$$ \frac{\rho_0}{g} \mathbb{D}_t b = (w^* dt + (\boldsymbol{\sigma}_t d\mathbf{B}_t)_z) \frac{\partial \rho_1}{\partial z} - \frac{1}{2} \nabla \cdot \left( \mathbf{a}_{\bullet z} \frac{\partial \rho_1}{\partial z} \right) dt, \quad (35) $$

with

$$ \mathbf{a} = \begin{pmatrix} \mathbf{a}_{HH^T} & \mathbf{a}_{Hz} \\ \mathbf{a}_{zH^T} & a_{zz} \end{pmatrix}, \quad (36) $$

and $\mathbf{a}_{zH^T} = \mathbf{a}_{Hz}^T$, for $H = (x \ y)^T$.

## Thermodynamic Effects

Equation (35) is part of the stochastic version of what is often referred to as the simple Boussinesq equations. In the ocean, thermodynamic effects can be important, and we propose to incorporate these effects by combining the buoyancy and the energy equation, following the steps of [22]. Assuming a linear equation of state for sea water, we have
# 5 Boussinesq-Hydrostatic Approximation

In this section, starting from the stochastic compressible Navier-Stokes equations, we perform the Boussinesq approximation by considering small density fluctuations. These fluctuations are neglected, when they are not multiplied by gravity $g$, which leads to the classical definition of the buoyancy. We perform as well the hydrostatic approximation through the classical aspect ratio scaling $D = H/L_{ref} \ll 1$, with $H$ the water depth. For simplicity, we do not consider a rotating frame. Coriolis correction could be straightforwardly considered as in [21]. The vertical coordinate $z \in [-H, \eta]$ is bounded by the bottom and the free surface.

## Density

The density is decomposed through the following asymptotic expansion

$$ \rho = \rho_0 + \epsilon\rho_1(z) + \epsilon\rho_2(x, y, z, t) + o(\epsilon), \quad (33) $$

with $\rho_1(z)$ the time-averaged stratification term, and $\epsilon \ll 1$ and we do not need to assume that $\rho_1 > \rho_2$. We obtain hence

$$ \nabla \cdot \mathbf{u}^* = 0 \quad ; \quad \nabla \cdot \boldsymbol{\sigma}_t \, d\mathbf{B}_t = 0 \quad ; \quad \mathbb{D}_t(\rho_1 + \rho_2) = 0. \quad (34) $$

The drift velocity and the noise are divergence free. Density perturbations undergo a stochastic transport by the flow. We remark that since $\nabla \cdot \boldsymbol{\sigma}_t d\mathbf{B}_t = 0$, then the transport operator $\mathbb{D}_t(\cdot)$ can be directly used.

The terms of order $\epsilon$ of Eq. (34) can be expressed in terms of buoyancy $b = -\epsilon g \rho_2 / \rho_0$:

$$ \frac{\rho_0}{g} \mathbb{D}_t b = (w^* dt + (\boldsymbol{\sigma}_t d\mathbf{B}_t)_z) \frac{\partial \rho_1}{\partial z} - \frac{1}{2} \nabla \cdot \left( \mathbf{a}_{\bullet z} \frac{\partial \rho_1}{\partial z} \right) dt, \quad (35) $$

with

$$ \mathbf{a} = \begin{pmatrix} \mathbf{a}_{HH^T} & \mathbf{a}_{Hz} \\ \mathbf{a}_{zH^T} & a_{zz} \end{pmatrix}, \quad (36) $$

and $\mathbf{a}_{zH^T} = \mathbf{a}_{Hz}^T$, for $H = (x \ y)^T$.

## Thermodynamic Effects

Equation (35) is part of the stochastic version of what is often referred to as the simple Boussinesq equations. In the ocean, thermodynamic effects can be important, and we propose to incorporate these effects by combining the buoyancy and the energy equation, following the steps of [22]. Assuming a linear equation of state for sea water, we have
$$ \rho = \rho_0(1 - \beta_T(T - T_0) + \beta_p p), \quad (37) $$

with $\beta_p = 1/\rho_0 c^2$ and $\beta_T = 1/\rho_0 \frac{\partial \rho}{\partial T}$ the coefficients of the Taylor expansion. For sake of simplicity, we do not take into account salinity effects, and we apply the stochastic transport operator to Eq. (37). We obtain

$$ \begin{aligned} \mathbb{D}_t \rho &= -\rho_0 \beta_T \mathbb{D}_t T + \frac{1}{c^2} \mathbb{D}_t p \\ \mathbb{D}_t \left( \rho - \frac{1}{c^2} p \right) &= -\frac{\beta_T}{\gamma} (dW + dQ). \end{aligned} \quad (38) $$

With no viscosity, divergence-free velocity and neglecting the quadratic variations (with the same argument as in (31)) together with the hydrostatic assumption on the leading term $p_0$, we can assume that the main dilatations are caused by radiative effects to which the potential buoyancy $b_\phi$ is directly sensitive:

$$ b_{\phi} \triangleq -\frac{g}{\rho_0} \left( \delta\rho + \frac{\rho_0 g z}{c^2} \right) = b_{st} + b - g \frac{z}{H_p}, \quad (39) $$

with $H_p = c^2/g$ and $b_{st} = -\epsilon g \rho_1 / \rho_0$. Upon applying the transport operator (with forcing), the following evolution equation of the potential buoyancy is obtained

$$ \mathbb{D}_t b_\phi = \frac{g \beta_T}{\gamma \rho_0} (dQ + dW_d), \quad (40) $$

where $dQ = dQ_{rad} + \frac{1}{RePr} \nabla \cdot (\nabla T) dt$, with $dQ_{rad}$ the radiative heat fluxes, and

$$ dW_d = (u^* - u) dt + \sigma_t d\mathbf{B}_t \cdot (-\nabla p + \rho g) - (u^* - u) \cdot \nabla p_t^\sigma $$

the drift works. The drift work on the vertical velocity component can be interpreted (with a linear equation of state) as an alternative to the so-called eddy diffusivity mass flux (EDMF) scheme proposed recently for atmospheric and oceanic penetrative convection parameterization (see for instance [19, 10] and references therein). Indeed, in EDMF, the subgrid stress in the transport equation of temperature is modelled as a mass flux induced by a given number of plumes (corresponding here possibly to $\rho \sigma_t d\mathbf{B}_t$), multiplied by the difference of temperature between the plume and the ambient flow, which is here proportional to a buoyancy anomaly. Interestingly, the pressure work provides a natural non-local (horizontal and vertical) forcing term while the other term is a local upward/downward vertical statistical forcing. EDMF schemes are obtained by specifying the noise in terms of velocity fluctuations between the mean velocity and non-convective environment, upward plumes and downward plumes. Such an interpretation need to be tested with numerical simulations, and will be the focus of a future dedicated study.
$$ \rho = \rho_0(1 - \beta_T(T - T_0) + \beta_p p), \quad (37) $$

with $\beta_p = 1/\rho_0 c^2$ and $\beta_T = 1/\rho_0 \frac{\partial \rho}{\partial T}$ the coefficients of the Taylor expansion. For sake of simplicity, we do not take into account salinity effects, and we apply the stochastic transport operator to Eq. (37). We obtain

$$ \begin{aligned} \mathbb{D}_t \rho &= -\rho_0 \beta_T \mathbb{D}_t T + \frac{1}{c^2} \mathbb{D}_t p \\ \mathbb{D}_t \left( \rho - \frac{1}{c^2} p \right) &= -\frac{\beta_T}{\gamma} (dW + dQ). \end{aligned} \quad (38) $$

With no viscosity, divergence-free velocity and neglecting the quadratic variations (with the same argument as in (31)) together with the hydrostatic assumption on the leading term $p_0$, we can assume that the main dilatations are caused by radiative effects to which the potential buoyancy $b_\phi$ is directly sensitive:

$$ b_{\phi} \triangleq -\frac{g}{\rho_0} \left( \delta\rho + \frac{\rho_0 g z}{c^2} \right) = b_{st} + b - g \frac{z}{H_p}, \quad (39) $$

with $H_p = c^2/g$ and $b_{st} = -\epsilon g \rho_1 / \rho_0$. Upon applying the transport operator (with forcing), the following evolution equation of the potential buoyancy is obtained

$$ \mathbb{D}_t b_\phi = \frac{g \beta_T}{\gamma \rho_0} (dQ + dW_d), \quad (40) $$

where $dQ = dQ_{rad} + \frac{1}{RePr} \nabla \cdot (\nabla T) dt$, with $dQ_{rad}$ the radiative heat fluxes, and

$$ dW_d = (u^* - u) dt + \sigma_t d\mathbf{B}_t \cdot (-\nabla p + \rho g) - (u^* - u) \cdot \nabla p_t^\sigma $$

the drift works. The drift work on the vertical velocity component can be interpreted (with a linear equation of state) as an alternative to the so-called eddy diffusivity mass flux (EDMF) scheme proposed recently for atmospheric and oceanic penetrative convection parameterization (see for instance [19, 10] and references therein). Indeed, in EDMF, the subgrid stress in the transport equation of temperature is modelled as a mass flux induced by a given number of plumes (corresponding here possibly to $\rho \sigma_t d\mathbf{B}_t$), multiplied by the difference of temperature between the plume and the ambient flow, which is here proportional to a buoyancy anomaly. Interestingly, the pressure work provides a natural non-local (horizontal and vertical) forcing term while the other term is a local upward/downward vertical statistical forcing. EDMF schemes are obtained by specifying the noise in terms of velocity fluctuations between the mean velocity and non-convective environment, upward plumes and downward plumes. Such an interpretation need to be tested with numerical simulations, and will be the focus of a future dedicated study.
By defining the buoyancy frequency

$$N^2(z) \triangleq \frac{\partial}{\partial z} \left( -\epsilon g \frac{\rho_1}{\rho_0} - g \frac{z}{H_p} \right) = -\epsilon \frac{g}{\rho_0} \frac{\partial \rho_1}{\partial z} - \frac{g^2}{c^2}, \quad (41)$$

stratification and radiative effects can be introduced explicitly on the buoyancy equation

$$\mathbb{D}_t b + (w^* dt + (\sigma_t d\mathbf{B}_t)_z) N^2 - \frac{1}{2} \nabla \cdot (\mathbf{a}_{\bullet z} N^2) dt = \frac{g \beta_T}{\gamma \rho_0} (dQ + dW_d). \quad (42)$$

## Momentum

Concerning the momentum equation, we neglect here the viscous terms. In this framework, g is assumed to be $O(1/\epsilon)$. We decompose as well the pressure field as follows

$$p = \underbrace{p_0(z)}_{O(1/\epsilon)} + p_1 + p_2 + O(\epsilon), \quad (43)$$

where $p_0$ and $p_1$ are in hydrostatic balance:

$$\frac{\partial p_0}{\partial z} = -g\rho_0 \quad \text{and} \quad \frac{\partial p_1}{\partial z} = -\epsilon g \rho_1. \quad (44)$$

The momentum equation (15) becomes

$$\begin{aligned} & \left( \rho_0 + \epsilon(\rho_1 + \rho_2) \right) \mathbb{D}_t u_i - \sum_k d_t \left\langle \int_0^t \rho \sigma_s d\mathbf{B}_s^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} \frac{\partial d p_s^\sigma}{\partial x_i} \right) \right\rangle \\ & = - \frac{\partial p_2}{\partial x_i} dt - \frac{\partial d p_t^\sigma}{\partial x_i} - \epsilon \rho_2 g \delta_{iz} dt. \end{aligned} \quad (45)$$

Similarly as in Sect. 4, the covariation term is $O(\epsilon)$. By neglecting $O(\epsilon)$ terms, we obtain

$$\mathbb{D}_t u_i = -\frac{1}{\rho_0} \frac{\partial p_2}{\partial x_i} dt - \frac{1}{\rho_0} \frac{\partial d p_t^\sigma}{\partial x_i} - \underbrace{-\epsilon \frac{\rho_2}{\rho_0} g \delta_{iz}}_{b} dt. \quad (46)$$

Finally, $p_2$ and the random pressure are determined through a generalisation of the hydrostatic balance, accounting for a part of non-hydrostatic effects by balancing in the vertical momentum equation the vertical pressure gradient with buoyancy, stochastic diffusion, corrective drift and stochastic advection of $w$. We consider a regime where the hydrostatic approximation in the deterministic framework is only roughly valid (in other words at the limit of validity), such that a noise with a strong amplitude can break this assumption—or changing viewpoint, the regime is intermediate and we aim at modelling some weak non-hydrostatic effects through stochastic modelling. By scaling analysis (weak aspect ratio and noise with strong amplitude), $d_t w$ and $(u \cdot \nabla) w$ are neglected while terms associated with the noise
By defining the buoyancy frequency

$$N^2(z) \triangleq \frac{\partial}{\partial z} \left( -\epsilon g \frac{\rho_1}{\rho_0} - g \frac{z}{H_p} \right) = -\epsilon \frac{g}{\rho_0} \frac{\partial \rho_1}{\partial z} - \frac{g^2}{c^2}, \quad (41)$$

stratification and radiative effects can be introduced explicitly on the buoyancy equation

$$\mathbb{D}_t b + (w^* dt + (\sigma_t d\mathbf{B}_t)_z) N^2 - \frac{1}{2} \nabla \cdot (\mathbf{a}_{\bullet z} N^2) dt = \frac{g \beta_T}{\gamma \rho_0} (dQ + dW_d). \quad (42)$$

## Momentum

Concerning the momentum equation, we neglect here the viscous terms. In this framework, g is assumed to be $O(1/\epsilon)$. We decompose as well the pressure field as follows

$$p = \underbrace{p_0(z)}_{O(1/\epsilon)} + p_1 + p_2 + O(\epsilon), \quad (43)$$

where $p_0$ and $p_1$ are in hydrostatic balance:

$$\frac{\partial p_0}{\partial z} = -g\rho_0 \quad \text{and} \quad \frac{\partial p_1}{\partial z} = -\epsilon g \rho_1. \quad (44)$$

The momentum equation (15) becomes

$$\begin{aligned} & \left( \rho_0 + \epsilon(\rho_1 + \rho_2) \right) \mathbb{D}_t u_i - \sum_k d_t \left\langle \int_0^t \rho \sigma_s d\mathbf{B}_s^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} \frac{\partial d p_s^\sigma}{\partial x_i} \right) \right\rangle \\ & = - \frac{\partial p_2}{\partial x_i} dt - \frac{\partial d p_t^\sigma}{\partial x_i} - \epsilon \rho_2 g \delta_{iz} dt. \end{aligned} \quad (45)$$

Similarly as in Sect. 4, the covariation term is $O(\epsilon)$. By neglecting $O(\epsilon)$ terms, we obtain

$$\mathbb{D}_t u_i = - \frac{1}{\rho_0} \frac{\partial p_2}{\partial x_i} dt - \frac{1}{\rho_0} \frac{\partial d p_t^\sigma}{\partial x_i} \underbrace{-\epsilon \frac{\rho_2}{\rho_0} g \delta_{iz}}_{b} dt. \quad (46)$$

Finally, $p_2$ and the random pressure are determined through a generalisation of the hydrostatic balance, accounting for a part of non-hydrostatic effects by balancing in the vertical momentum equation the vertical pressure gradient with buoyancy, stochastic diffusion, corrective drift and stochastic advection of $w$. We consider a regime where the hydrostatic approximation in the deterministic framework is only roughly valid (in other words at the limit of validity), such that a noise with a strong amplitude can break this assumption—or changing viewpoint, the regime is intermediate and we aim at modelling some weak non-hydrostatic effects through stochastic modelling. By scaling analysis (weak aspect ratio and noise with strong amplitude), $d_t w$ and $(u \cdot \nabla) w$ are neglected while terms associated with the noise
are kept. Indeed, denoting $L^\sigma$ the scale amplitude of $\sigma_t d\mathbf{B}_t$, and $\tau$ the decorrelation time, the advection of $\sigma_t d\mathbf{B}_t$ cannot be neglected if<fnref n="1" /> $L^\sigma/L_{\text{ref}} \sim 1/(Fr D)^2$ and stochastic diffusion and drift velocity are important if $(L^\sigma/L_{\text{ref}})^2\tau/T_{\text{ref}} \sim 1/(Fr D)^2$, with the Froude number $Fr = u_{\text{ref}}/(NH)$. Since $d_t w$ is neglected, martingale and time-differentiable terms can then be safely separated, such that the remaining pressure term can be determined by vertical integration through a scheme similar to the one applied in the classical hydrostatic regime:

$$ 
\begin{align} 
p_2 &= \rho_0 \int_z^\eta \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) - b \right) dz \tag{47} \\ 
d p_t^\sigma &= -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz. \nonumber 
\end{align} 
$$

Here, we have neglected $d_t w$, but random vertical transport could generate some random vertical acceleration instead of only random pressure fluctuations. An intermediary assumption could be to consider that the time-differentiable part of $d_t w$ is negligible (classical hydrostatic balance), but that its martingale part is not. It could be obtained by diagnosing the vertical velocity time increment and thus bringing and additional correction to $d p_t^\sigma$ in Eq. (47).

## Summary

By collecting the Eqs. (42), (46), and (47), we obtain the following stochastic Boussinesq system with thermodynamic forcing

$$ 
\left\{ 
\begin{aligned} 
\mathbb{D}_t u_i &= -\frac{1}{\rho_0} \frac{\partial p_2}{\partial x_i} dt - \frac{1}{\rho_0} \frac{\partial d p_t^\sigma}{\partial x_i} \quad \text{for } i = \{u, v\} \\ 
w &= \frac{1}{2} (\nabla \cdot \mathbf{a})_z - \int_{-H}^z \left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right) dz \\ 
\nabla \cdot (\sigma_t d\mathbf{B}_t) &= 0 \\ 
p_2 &= \rho_0 \int_z^\eta \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) - b \right) dz \\ 
d p_t^\sigma &= -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz \\ 
\mathbb{D}_t b + (w^* dt + (\sigma_t d\mathbf{B}_t)_z) N^2 - \frac{1}{2} \nabla \cdot (a_\bullet z N^2) \, dt &= \frac{g\beta_T}{\gamma\rho_0} (dQ + dW_d) \\ 
dQ &= dQ_{\text{rad}} + \frac{1}{RePr} \nabla \cdot (\nabla T) \, dt \\ 
dW_d &= \left( \frac{1}{2} \nabla \cdot \mathbf{a} \, dt - \sigma_t d\mathbf{B}_t \right) \cdot (\nabla p_2 + \rho_0 b e_z) + \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla d p_t^\sigma. 
\end{aligned} 
\right. 
\tag{48} 
$$

<footnote n="1">¹ If the frame rotation is taken into account, the ratio $R_o/Bu$ Rossby over Burger is additionally involved, but does not change the existence of an intermediary regime.</footnote>
are kept. Indeed, denoting $L^\sigma$ the scale amplitude of $\sigma_t d\mathbf{B}_t$, and $\tau$ the decorrelation time, the advection of $\sigma_t d\mathbf{B}_t$ cannot be neglected if<fnref n="1" /> $L^\sigma/L_{\text{ref}} \sim 1/(Fr D)^2$ and stochastic diffusion and drift velocity are important if $(L^\sigma/L_{\text{ref}})^2\tau/T_{\text{ref}} \sim 1/(Fr D)^2$, with the Froude number $Fr = u_{\text{ref}}/(NH)$. Since $d_t w$ is neglected, martingale and time-differentiable terms can then be safely separated, such that the remaining pressure term can be determined by vertical integration through a scheme similar to the one applied in the classical hydrostatic regime:

$$ 
\begin{align} 
p_2 &= \rho_0 \int_z^\eta \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) - b \right) dz \tag{47} \\ 
d p_t^\sigma &= -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz. \nonumber
\end{align} 
$$

Here, we have neglected $d_t w$, but random vertical transport could generate some random vertical acceleration instead of only random pressure fluctuations. An intermediary assumption could be to consider that the time-differentiable part of $d_t w$ is negligible (classical hydrostatic balance), but that its martingale part is not. It could be obtained by diagnosing the vertical velocity time increment and thus bringing and additional correction to $d p_t^\sigma$ in Eq. (47).

## Summary

By collecting the Eqs. (42), (46), and (47), we obtain the following stochastic Boussinesq system with thermodynamic forcing

$$ 
\left\{ 
\begin{aligned}
\mathbb{D}_t u_i &= -\frac{1}{\rho_0} \frac{\partial p_2}{\partial x_i} dt - \frac{1}{\rho_0} \frac{\partial d p_t^\sigma}{\partial x_i} \quad \text{for } i = \{u, v\} \\
w &= \frac{1}{2} (\nabla \cdot \mathbf{a})_z - \int_{-H}^z \left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right) dz \\
\nabla \cdot (\sigma_t d\mathbf{B}_t) &= 0 \\
p_2 &= \rho_0 \int_z^\eta \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) - b \right) dz \\
d p_t^\sigma &= -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz \\
\mathbb{D}_t b + (w^* dt + (\sigma_t d\mathbf{B}_t)_z) N^2 - \frac{1}{2} \nabla \cdot (a_\bullet z N^2) \, dt &= \frac{g\beta_T}{\gamma\rho_0} (dQ + dW_d) \\
dQ &= dQ_{\text{rad}} + \frac{1}{RePr} \nabla \cdot (\nabla T) \, dt \\
dW_d &= \left( \frac{1}{2} \nabla \cdot \mathbf{a} \, dt - \sigma_t d\mathbf{B}_t \right) \cdot (\nabla p_2 + \rho_0 b e_z) + \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla d p_t^\sigma.
\end{aligned} 
\right. 
\tag{48} 
$$

<footnote n="1">¹ If the frame rotation is taken into account, the ratio $R_o/Bu$ Rossby over Burger is additionally involved, but does not change the existence of an intermediary regime.</footnote>
In system (48), the pressure is obtained through a relaxed hydrostatic balance, and the vertical velocity is deduced kinematically from the divergence-free condition of the drift velocity. Neglecting the thermodynamic effects, together with a strong hydrostatic balance assumption (weak to moderate noise regime), we recover the simple Boussinesq system presented in [16], without the Coriolis correction. Obviously, this latter could be added without any major difficulty.

In some applications, a more accurate evaluation of the buoyancy is required, and it can be obtained through an equation of state $\rho_{\text{BQ}}(T, p)$ (salinity is not taken into account here and left for future works) associated with a transport equation of temperature (and salinity when considered). Under the aforementioned assumptions, the transport equation of temperature (23) is simplified, and the full system can be written

$$ 
\left\{
\begin{aligned}
\mathbb{D}_t u_i &= -\frac{1}{\rho_0} \frac{\partial p_2}{\partial x_i} dt - \frac{1}{\rho_0} \frac{\partial d p_t^\sigma}{\partial x_i} \quad \text{for } i = \{u, v\} \\
w &= \frac{1}{2} (\nabla \cdot \mathbf{a})_z - \int_{-H}^{z} \left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right) dz \\
\nabla \cdot (\boldsymbol{\sigma}_t d\mathbf{B}_t) &= 0 \\
p_2 &= \rho_0 \int_z^\eta \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) - b \right) dz \\
d p_t^\sigma &= -\rho_0 \int_z^\eta (\boldsymbol{\sigma}_t d\mathbf{B}_t \cdot \nabla) w \, dz \\
\frac{\rho}{\gamma} \mathbb{D}_t T &= \left( \frac{1}{2} \nabla \cdot \mathbf{a} \, dt - \boldsymbol{\sigma}_t d\mathbf{B}_t \right) \cdot (\nabla p_2 + \rho_0 b e_z) + \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla d p_t^\sigma \\
&\quad + \frac{1}{RePr} \nabla \cdot (\nabla T) \, dt + dQ_{\text{rad}} \\
b &= -\frac{g}{\rho_0} \rho_{\text{BQ}}(T, -\rho_0 gz).
\end{aligned}
\right.
\quad (49) $$

Usually, a simple stochastic advection-diffusion equation is considered for the transport of temperature, but in the system (49), we can point out that the drift works remain. These source/sink terms in the temperature evolution equation is one of the principal outcome of this study. As outlined this additional terms for parameterising discrepancies to hydrostatic physics and primitive equations. In the next section, we explore systems at finer resolution.
In system (48), the pressure is obtained through a relaxed hydrostatic balance, and the vertical velocity is deduced kinematically from the divergence-free condition of the drift velocity. Neglecting the thermodynamic effects, together with a strong hydrostatic balance assumption (weak to moderate noise regime), we recover the simple Boussinesq system presented in [16], without the Coriolis correction. Obviously, this latter could be added without any major difficulty.

In some applications, a more accurate evaluation of the buoyancy is required, and it can be obtained through an equation of state $\rho_{\text{BQ}}(T, p)$ (salinity is not taken into account here and left for future works) associated with a transport equation of temperature (and salinity when considered). Under the aforementioned assumptions, the transport equation of temperature (23) is simplified, and the full system can be written

$$ 
\left\{
\begin{aligned}
\mathbb{D}_t u_i &= -\frac{1}{\rho_0} \frac{\partial p_2}{\partial x_i} dt - \frac{1}{\rho_0} \frac{\partial d p_t^\sigma}{\partial x_i} \quad \text{for } i = \{u, v\} \\
w &= \frac{1}{2} (\nabla \cdot \mathbf{a})_z - \int_{-H}^{z} \left( \frac{\partial u^*}{\partial x} + \frac{\partial v^*}{\partial y} \right) dz \\
\nabla \cdot (\boldsymbol{\sigma}_t d\mathbf{B}_t) &= 0 \\
p_2 &= \rho_0 \int_z^\eta \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) - b \right) dz \\
d p_t^\sigma &= -\rho_0 \int_z^\eta (\boldsymbol{\sigma}_t d\mathbf{B}_t \cdot \nabla) w \, dz \\
\frac{\rho}{\gamma} \mathbb{D}_t T &= \left( \frac{1}{2} \nabla \cdot \mathbf{a} \, dt - \boldsymbol{\sigma}_t d\mathbf{B}_t \right) \cdot (\nabla p_2 + \rho_0 b e_z) + \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla d p_t^\sigma \\
&\quad + \frac{1}{RePr} \nabla \cdot (\nabla T) \, dt + dQ_{\text{rad}} \\
b &= -\frac{g}{\rho_0} \rho_{\text{BQ}}(T, -\rho_0 g z).
\end{aligned}
\right.
\quad (49) $$

Usually, a simple stochastic advection-diffusion equation is considered for the transport of temperature, but in the system (49), we can point out that the drift works remain. These source/sink terms in the temperature evolution equation is one of the principal outcome of this study. As outlined this additional terms for parameterising discrepancies to hydrostatic physics and primitive equations. In the next section, we explore systems at finer resolution.
# 6 Extension to Non-Boussinesq

The aim of this section is to propose a formulation to relax the Boussinesq assumption in the LU stochastic framework while avoiding the resolution of a 3D Poisson equation. We consider now an intermediate model between the fully non-Boussinesq non-hydrostatic formulation and the system (49). For sea water, the equation of state is formulated in terms of $\rho_{BQ}(T, p)$ instead of $p(\rho, T)$ as in gas dynamics. To take this aspect into consideration, we follow [3] in order to obtain an explicit expression of the pressure. This is at the cost of resolving in time sound waves, or a pseudo-compressibility information propagating at the velocity $c$. The density is decomposed as

$$ \rho = \rho_{BQ}(T, p) + \underbrace{\frac{\partial \rho}{\partial p} \delta p}_{\delta \rho} + \mathcal{O}(\delta p^2), \quad (50) $$

with $\rho_{BQ}(T, -\rho_0 g z)$ the Boussinesq density determined by the equation of state under an hydrostatic balance condition. The deviation to this density is then assumed to be ensue from an isentropic transformation, i.e. of acoustic nature. The term $\frac{\partial \rho}{\partial p} = \frac{1}{c^2}$, is then directly related to the sound speed (or more precisely to the fastest wave considered in the model). We determine now a transport equation for $\delta \rho$.

## Continuity

We start from the continuity equation of the stochastic compressible Navier-Stokes equations

$$ d_t \rho + \nabla \cdot \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t d\mathbf{B}_t \right) \rho = \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \rho) dt. \quad (51) $$

A transport equation for $\delta \rho$ can be deduced as

$$ d_t (\delta \rho) = -d_t (\rho_{BQ}) - \nabla \cdot \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t d\mathbf{B}_t \right) (\rho_{BQ} + \delta \rho) + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \rho_{BQ} + \delta \rho) dt. \quad (52) $$

At the scales considered here we assume that the unresolved contribution is of hydrodynamic nature associated to a divergence free noise $\nabla \cdot (\sigma_t d\mathbf{B}_t) = 0$.

## Momentum

The pressure can be decomposed as well as

$$ p = p_{\text{atm}} + \int_z^\eta \rho_{\text{BQ}}(z') g \, dz' + p_{\text{NH}} + c^2 \delta \rho, \quad (53) $$
# 6 Extension to Non-Boussinesq

The aim of this section is to propose a formulation to relax the Boussinesq assumption in the LU stochastic framework while avoiding the resolution of a 3D Poisson equation. We consider now an intermediate model between the fully non-Boussinesq non-hydrostatic formulation and the system (49). For sea water, the equation of state is formulated in terms of $\rho_{BQ}(T, p)$ instead of $p(\rho, T)$ as in gas dynamics. To take this aspect into consideration, we follow [3] in order to obtain an explicit expression of the pressure. This is at the cost of resolving in time sound waves, or a pseudo-compressibility information propagating at the velocity $c$. The density is decomposed as

$$ \rho = \rho_{BQ}(T, p) + \underbrace{\frac{\partial \rho}{\partial p} \delta p}_{\delta \rho} + \mathcal{O}(\delta p^2), \quad (50) $$

with $\rho_{BQ}(T, -\rho_0 g z)$ the Boussinesq density determined by the equation of state under an hydrostatic balance condition. The deviation to this density is then assumed to be ensue from an isentropic transformation, i.e. of acoustic nature. The term $\frac{\partial \rho}{\partial p} = \frac{1}{c^2}$, is then directly related to the sound speed (or more precisely to the fastest wave considered in the model). We determine now a transport equation for $\delta \rho$.

## Continuity

We start from the continuity equation of the stochastic compressible Navier-Stokes equations

$$ d_t \rho + \nabla \cdot \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t d\mathbf{B}_t \right) \rho = \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \rho) dt. \quad (51) $$

A transport equation for $\delta \rho$ can be deduced as

$$ d_t(\delta\rho) = -d_t(\rho_{BQ}) - \nabla \cdot \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t d\mathbf{B}_t \right) (\rho_{BQ} + \delta\rho) + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \rho_{BQ} + \delta\rho) dt. \quad (52) $$

At the scales considered here we assume that the unresolved contribution is of hydrodynamic nature associated to a divergence free noise $\nabla \cdot (\sigma_t d\mathbf{B}_t) = 0$.

## Momentum

The pressure can be decomposed as well as

$$ p = p_{\text{atm}} + \int_{z}^{\eta} \rho_{\text{BQ}}(z') g \, dz' + p_{\text{NH}} + c^2 \delta \rho, \quad (53) $$
where $p_{atm}$ is the atmospheric pressure, which will be neglected later for simplicity. The second term is the hydrostatic pressure associated with the Boussinesq density. The third term, $p_{NH}$, is associated with Boussinesq non-hydrostatic effects balancing vertical advection, and finally $c^2\delta\rho$ corresponds to a non-Boussinesq component of acoustic nature.

For the martingale random pressure, two components are considered: a Boussinesq non-hydrostatic term as in Eq. (48), and a non-Boussinesq component of acoustic nature as in Eq. (27)

$$d p_t^\sigma = -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz' - \tau c^2 (\sigma_t d\mathbf{B}_t \cdot \nabla) \rho. \quad (54)$$

Neglecting the viscous terms, we obtain for the momentum equation

$$\rho \mathbb{D}_t u_i - \sum_k \mathbb{D}_t \left\langle \int_0^t \rho (\sigma_s d\mathbf{B}_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} \frac{\partial d p_s^\sigma}{\partial x_i} \right) \right\rangle = -\frac{\partial p}{\partial x_i} \, dt - \frac{\partial d p_t^\sigma}{\partial x_i} + \rho \mathbf{g}. \quad (55)$$

Assuming that $\rho \mathbb{D}_t u_i \approx \rho_{BQ} \mathbb{D}_t u_i$, and following the same arguments as in Sect. 4 to neglect the quadratic variation term, one finally get

$$\rho_{BQ} \mathbb{D}_t u = -\frac{\partial p}{\partial x} dt - \frac{\partial d p_t^\sigma}{\partial x}$$

$$\rho_{BQ} \mathbb{D}_t v = -\frac{\partial p}{\partial y} dt - \frac{\partial d p_t^\sigma}{\partial y}$$

$$\rho_{BQ} \mathbb{D}_t w = -\frac{\partial p}{\partial z} dt - \frac{\partial d p_t^\sigma}{\partial z} + (\rho_{BQ} + \delta\rho) \mathbf{g} dt$$

$$p = \int_z^\eta \left( (\rho_{BQ}(z') + \delta\rho)g + \rho_{BQ} \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) \right) \right) dz' + c^2 \delta\rho$$

$$d p_t^\sigma = -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz' - \tau c^2 (\sigma_t d\mathbf{B}_t \cdot \nabla) \rho. \quad (56)$$

The system (52)–(56) can be solved explicitly and does not require the expensive resolution of a 3D Poisson equation. Although system (49) proposes a deviation to the hydrostatic hypothesis through the martingale random pressure, the system (56) considers a non-hydrostatic model that fully accounts for stochastic vertical accelerations while relaxing the effect of fast waves truncation through the martingale pressure term. This system remains restricted by a CFL condition depending on the propagation speed of pseudo-compressibility informations. We believe this modelling strategy opens some new research directions on the role of unresolved small scales on non-hydrostatic and non-Boussinesq effects in oceanic flows.
where $p_{atm}$ is the atmospheric pressure, which will be neglected later for simplicity. The second term is the hydrostatic pressure associated with the Boussinesq density. The third term, $p_{NH}$, is associated with Boussinesq non-hydrostatic effects balancing vertical advection, and finally $c^2\delta\rho$ corresponds to a non-Boussinesq component of acoustic nature.

For the martingale random pressure, two components are considered: a Boussinesq non-hydrostatic term as in Eq. (48), and a non-Boussinesq component of acoustic nature as in Eq. (27)

$$d p_t^\sigma = -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz' - \tau c^2 (\sigma_t d\mathbf{B}_t \cdot \nabla) \rho. \quad (54)$$

Neglecting the viscous terms, we obtain for the momentum equation

$$\rho \mathbb{D}_t u_i - \sum_k \mathbb{D}_t \left\langle \int_0^t \rho (\sigma_s d\mathbf{B}_s)^k, \int_0^t \frac{\partial}{\partial x_k} \left( \frac{1}{\rho} \frac{\partial d p_s^\sigma}{\partial x_i} \right) \right\rangle = -\frac{\partial p}{\partial x_i} \, dt - \frac{\partial d p_t^\sigma}{\partial x_i} + \rho \mathbf{g}. \quad (55)$$

Assuming that $\rho \mathbb{D}_t u_i \approx \rho_{BQ} \mathbb{D}_t u_i$, and following the same arguments as in Sect. 4 to neglect the quadratic variation term, one finally get

$$\rho_{BQ} \mathbb{D}_t u = -\frac{\partial p}{\partial x} dt - \frac{\partial d p_t^\sigma}{\partial x}$$

$$\rho_{BQ} \mathbb{D}_t v = -\frac{\partial p}{\partial y} dt - \frac{\partial d p_t^\sigma}{\partial y}$$

$$\rho_{BQ} \mathbb{D}_t w = -\frac{\partial p}{\partial z} dt - \frac{\partial d p_t^\sigma}{\partial z} + (\rho_{BQ} + \delta\rho) \mathbf{g} dt$$

$$p = \int_z^\eta \left( (\rho_{BQ}(z') + \delta\rho)g + \rho_{BQ} \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \cdot \nabla \right) w + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) \right) \right) dz' + c^2 \delta\rho$$

$$d p_t^\sigma = -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \, dz' - \tau c^2 (\sigma_t d\mathbf{B}_t \cdot \nabla) \rho. \quad (56)$$

The system (52)–(56) can be solved explicitly and does not require the expensive resolution of a 3D Poisson equation. Although system (49) proposes a deviation to the hydrostatic hypothesis through the martingale random pressure, the system (56) considers a non-hydrostatic model that fully accounts for stochastic vertical accelerations while relaxing the effect of fast waves truncation through the martingale pressure term. This system remains restricted by a CFL condition depending on the propagation speed of pseudo-compressibility informations. We believe this modelling strategy opens some new research directions on the role of unresolved small scales on non-hydrostatic and non-Boussinesq effects in oceanic flows.
# 7 Conclusion

This paper proposes a stochastic representation under location uncertainty of the compressible Navier-Stokes equations. It as been obtained from conservation of density, momentum and total energy, undergoing a stochastic transport. The structure of equations remains similar to the compressible deterministic case. Nevertheless, because of the specificities related to stochastic transport, we have identified additional terms such as work induced by the alignment between the time-differentiable pressure gradient and the drift velocity. This small scale induced work is alike the baropycnal work known in compressible large eddy simulations and includes also terms reminiscent to mass flux parameterisation of atmospheric and oceanic penetrative convection phenomenon. These terms are obtained by the mean of a rigorous derivation from the conservation laws coupled with stochastic calculus rules associated to stochastic transport, instead of phenomenological arguments.

We have verified that applying low-Mach and Boussinesq approximations on the stochastic compressible system enabled us to recover the known incompressible and simple Boussinesq stochastic systems respectively. The general set of stochastic compressible equations allowed us to incorporate thermodynamic effects on the Boussinesq system. Finally, this formulation has lead us to propose a way to relax the Boussinesq and hydrostatic assumptions. This study opens some new research directions to exploit the potential of stochastic modelling for the numerical simulations of oceanic flows we will exploit in future works.

**Acknowledgments** The authors acknowledge the support of the ERC EU project 856408-STUOD and warmly thank Valentin Resseguier for very fruitful and stimulating discussions on this work.

# Appendix A: Stochastic Reynolds Transport Theorem from Stratonovich to Itō

The aim of this section is to rewrite the stochastic Reynolds transport theorem with a Stratonovich convention, and verify that we get the same relation as the ones obtained in Sect. 2 in the Itō setting. To that end, we follow the steps of [16, Appendix D] but we use the Stratonovich convention. Then, we pass from the Stratonovich to the Itō form. These forms are equivalent for regular enough processes, thus allowing us to verify the consistency of the Eq. (9).

To pass from Itō to Stratonovich, we use the following relation

$$X_t \circ dY_t = X_t dY_t + \frac{1}{2} d_t \left\langle \int_0^t dX_s, \int_0^t dY_s \right\rangle. \quad (57)$$

We define the characteristic function $\phi(x, t)$ transported by the flow, such that
# 7 Conclusion

This paper proposes a stochastic representation under location uncertainty of the compressible Navier-Stokes equations. It as been obtained from conservation of density, momentum and total energy, undergoing a stochastic transport. The structure of equations remains similar to the compressible deterministic case. Nevertheless, because of the specificities related to stochastic transport, we have identified additional terms such as work induced by the alignment between the time-differentiable pressure gradient and the drift velocity. This small scale induced work is alike the baropycnal work known in compressible large eddy simulations and includes also terms reminiscent to mass flux parameterisation of atmospheric and oceanic penetrative convection phenomenon. These terms are obtained by the mean of a rigorous derivation from the conservation laws coupled with stochastic calculus rules associated to stochastic transport, instead of phenomenological arguments.

We have verified that applying low-Mach and Boussinesq approximations on the stochastic compressible system enabled us to recover the known incompressible and simple Boussinesq stochastic systems respectively. The general set of stochastic compressible equations allowed us to incorporate thermodynamic effects on the Boussinesq system. Finally, this formulation has lead us to propose a way to relax the Boussinesq and hydrostatic assumptions. This study opens some new research directions to exploit the potential of stochastic modelling for the numerical simulations of oceanic flows we will exploit in future works.

**Acknowledgments** The authors acknowledge the support of the ERC EU project 856408-STUOD and warmly thank Valentin Resseguier for very fruitful and stimulating discussions on this work.

# Appendix A: Stochastic Reynolds Transport Theorem from Stratonovich to Itō

The aim of this section is to rewrite the stochastic Reynolds transport theorem with a Stratonovich convention, and verify that we get the same relation as the ones obtained in Sect. 2 in the Itō setting. To that end, we follow the steps of [16, Appendix D] but we use the Stratonovich convention. Then, we pass from the Stratonovich to the Itō form. These forms are equivalent for regular enough processes, thus allowing us to verify the consistency of the Eq. (9).

To pass from Itō to Stratonovich, we use the following relation

$$X_t \circ dY_t = X_t dY_t + \frac{1}{2} d_t \left\langle \int_0^t dX_s, \int_0^t dY_s \right\rangle. \quad (57)$$

We define the characteristic function $\phi(x, t)$ transported by the flow, such that
$$ \phi(\mathbf{X}_t(\mathbf{x}_0)) = g(\mathbf{x}_0), \qquad (58) $$

with a compact spatial support $V(t)$ of non-zero values that does not include points on the domain boundary. We can then write

$$ \begin{aligned} d \int_{V(t)} (q\phi)(\mathbf{x}, t) \, d\mathbf{x} &= d \int_{\Omega} (q\phi)(\mathbf{x}, t) \, d\mathbf{x} \tag{59} \\ &= \int_{\Omega} d_t \circ q \, \phi + q \, d_t \circ \phi \, d\mathbf{x}. \end{aligned} $$

Since $\phi$ is transported, we have

$$ \begin{aligned} d\phi(\mathbf{X}_t, t) &= d_t \circ \phi + \nabla\phi \cdot d\mathbf{X}_t = 0 \tag{60} \\ d_t \circ \phi + \left(u - \frac{1}{2}\nabla \cdot \mathbf{a} + \frac{1}{2}\sigma_t(\nabla \cdot \mathbf{\sigma}_t)\right) \cdot \nabla\phi + (\nabla\phi \cdot \mathbf{\sigma}_t) \circ d\mathbf{B}_t &= 0. \end{aligned} $$

We have then

$$ \begin{aligned} & d \int_{V(t)} (q\phi)(\mathbf{x}, t) \, d\mathbf{x} \\ &= \int_{\Omega} d_t \circ q \, \phi - q \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \sigma_t (\nabla \cdot \mathbf{\sigma}_t) \right) \cdot \nabla \phi + (\nabla \phi \cdot \mathbf{\sigma}_t) \circ d\mathbf{B}_t \right) \, d\mathbf{x} \\ &= \int_{\Omega} \left[ d_t \circ q + \nabla \cdot \left( q \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \sigma_t (\nabla \cdot \mathbf{\sigma}_t) \right) dt + \mathbf{\sigma}_t \circ d\mathbf{B}_t \right) \right) \right] \phi \, d\mathbf{x}. \end{aligned} \qquad (61) $$

We add now a force and obtain the SRTT in Stratonovich form:

$$ d_t \circ q + \nabla \cdot \left( q \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \sigma_t (\nabla \cdot \mathbf{\sigma}_t) \right) dt + \mathbf{\sigma}_t \circ d\mathbf{B}_t \right) \right) = Q_t dt + Q_\sigma \circ d\mathbf{B}_t. \qquad (62) $$

Let us now write this in Itō form:

$$ \begin{aligned} \frac{\partial}{\partial x_i} (q \sigma_{t,ij} \circ d\mathbf{B}_t^j) = & \frac{\partial}{\partial x_i} (q \sigma_{t,ij} d\mathbf{B}_t^j) \\ & + \frac{1}{2} d_t \left\langle \int_0^t d_t \left( \frac{\partial}{\partial x_i} (q \sigma_{s,ij}) \right), \int_0^t d\mathbf{B}_s^j \right\rangle. \end{aligned} \qquad (63) $$

Since $\sigma_t$ is time differentiable in the Eulerian grid, we have

$$ d_t \left\langle \int_0^t d_t \sigma_{s,ij}, \int_0^t d\mathbf{B}_s^j \right\rangle = 0. $$
$$ \phi(\mathbf{X}_t(\mathbf{x}_0)) = g(\mathbf{x}_0), \qquad (58) $$

with a compact spatial support $V(t)$ of non-zero values that does not include points on the domain boundary. We can then write

$$ \begin{aligned} d \int_{V(t)} (q\phi)(\mathbf{x}, t) \, d\mathbf{x} &= d \int_{\Omega} (q\phi)(\mathbf{x}, t) \, d\mathbf{x} \tag{59} \\ &= \int_{\Omega} d_t \circ q \, \phi + q \, d_t \circ \phi \, d\mathbf{x}. \end{aligned} $$

Since $\phi$ is transported, we have

$$ \begin{aligned} d\phi(\mathbf{X}_t, t) &= d_t \circ \phi + \nabla\phi \cdot d\mathbf{X}_t = 0 \tag{60} \\ d_t \circ \phi + \left(u - \frac{1}{2}\nabla \cdot \mathbf{a} + \frac{1}{2}\sigma_t(\nabla \cdot \mathbf{\sigma}_t)\right) \cdot \nabla\phi + (\nabla\phi \cdot \mathbf{\sigma}_t) \circ d\mathbf{B}_t &= 0. \end{aligned} $$

We have then

$$ \begin{aligned} & d \int_{V(t)} (q\phi)(\mathbf{x}, t) \, d\mathbf{x} \\ &= \int_{\Omega} d_t \circ q \, \phi - q \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \sigma_t (\nabla \cdot \mathbf{\sigma}_t) \right) \cdot \nabla \phi + (\nabla \phi \cdot \mathbf{\sigma}_t) \circ d\mathbf{B}_t \right) \, d\mathbf{x} \\ &= \int_{\Omega} \left[ d_t \circ q + \nabla \cdot \left( q \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \sigma_t (\nabla \cdot \mathbf{\sigma}_t) \right) dt + \mathbf{\sigma}_t \circ d\mathbf{B}_t \right) \right) \right] \phi \, d\mathbf{x}. \end{aligned} \qquad (61) $$

We add now a force and obtain the SRTT in Stratonovich form:

$$ d_t \circ q + \nabla \cdot \left( q \left( \left( u - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \sigma_t (\nabla \cdot \mathbf{\sigma}_t) \right) dt + \mathbf{\sigma}_t \circ d\mathbf{B}_t \right) \right) = Q_t dt + Q_\sigma \circ d\mathbf{B}_t. \qquad (62) $$

Let us now write this in Itō form:

$$ \begin{aligned} \frac{\partial}{\partial x_i} (q \sigma_{t,ij} \circ d\mathbf{B}_t^j) = & \frac{\partial}{\partial x_i} (q \sigma_{t,ij} d\mathbf{B}_t^j) \\ & + \frac{1}{2} d_t \left\langle \int_0^t d_t \left( \frac{\partial}{\partial x_i} (q \sigma_{s,ij}) \right), \int_0^t d\mathbf{B}_s^j \right\rangle. \end{aligned} \qquad (63) $$

Since $\sigma_t$ is time differentiable in the Eulerian grid, we have

$$ d_t \left\langle \int_0^t d_t \sigma_{s,ij}, \int_0^t d\mathbf{B}_s^j \right\rangle = 0. $$
Then,

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

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

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

We can then write

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

Assembling everything and dropping the space integral, we obtain
Then,

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

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

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

We can then write

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

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

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

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

which is exactly Eq. (9).

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

## Appendix B: Calculation Rules

### Distributivity of the Stochastic Transport Operator

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

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

we have then

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

or less formally:
$$ 
\begin{aligned}
& d_t \circ q + \nabla \cdot \left( q \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \boldsymbol{\sigma}_t (\nabla \cdot \boldsymbol{\sigma}_t) \right) dt + \boldsymbol{\sigma}_t \circ d\mathbf{B}_t \right) \right) \\
& - Q_t \, dt - Q_\sigma \circ d\mathbf{B}_t \\
&= d_t q + \nabla \cdot \left( q \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} + \frac{1}{2} \boldsymbol{\sigma}_t (\nabla \cdot \boldsymbol{\sigma}_t) \right) dt + \boldsymbol{\sigma}_t d\mathbf{B}_t \right) \right) \\
& \quad + \frac{1}{2} \left[ \nabla \cdot (\boldsymbol{\sigma}_t Q_\sigma) dt - \nabla \cdot (\mathbf{a} \nabla q) dt - \nabla \cdot (q \boldsymbol{\sigma}_t (\nabla \cdot \boldsymbol{\sigma}_t)) dt \right] \\
& \quad - Q_t \, dt - Q_\sigma d\mathbf{B}_t + \frac{1}{2} \nabla \cdot (\boldsymbol{\sigma}_t Q_\sigma) \, dt.
\end{aligned}
\quad (67) $$

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

$$ 
\begin{aligned}
& d_t q + \nabla \cdot \left( q \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \boldsymbol{\sigma}_t d\mathbf{B}_t \right) \right) + \nabla \cdot (\boldsymbol{\sigma}_t Q_\sigma) dt \\
&= \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla q) dt + Q_t dt + Q_\sigma d\mathbf{B}_t,
\end{aligned}
\quad (68) $$

which is exactly Eq. (9).

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

## Appendix B: Calculation Rules

### Distributivity of the Stochastic Transport Operator

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

$$ 
\begin{aligned}
\mathbb{D}_t f &= F_t \, dt + \mathbf{F}_\sigma \cdot d\mathbf{B}_t \\
\mathbb{D}_t g &= G_t \, dt + \mathbf{G}_\sigma \cdot d\mathbf{B}_t,
\end{aligned}
\quad (69) $$

we have then

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

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

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

Proof

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

Developing only the covariation term:

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

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

Proof

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

Developing only the covariation term:

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

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

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

## Work of Random Forces

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

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

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

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

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

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

by taking $\psi(x, t) = \phi(t)\mathbf{u}(x, t)$, we obtain
whose right hand side is written formally

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

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

## Work of Random Forces

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

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

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

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

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

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

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

We can identify

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

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

## Appendix C: Displacement of a Transported Control Surface

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

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

Hence, the normal displacement of the control surface is $dX_{d,t} \cdot n$, which involves the modified advection velocity. As a consequence, the modified advection velocity has to be considered for the definitions of elementary works based on surface integrals.
$$ 
\begin{aligned}
\int_0^T \phi(t) dW_{\sigma,u} &= \int_0^T \phi(t) \mathbf{u}(x,t) d\mathbf{F}_{\sigma} - \int_0^T d\left(\phi(t) \mathbf{u}(x,t) \cdot \int_0^t d\mathbf{F}_{\sigma}\right) \\
&= \int_0^T \phi(t) \mathbf{u}(x,t) d\mathbf{F}_{\sigma} - \phi(T) \mathbf{u}(x,T) \cdot \mathbf{F}_{\sigma}(x,T) \quad (78) \\
&= \int_0^T \phi(t) \mathbf{u}(x,t) d\mathbf{F}_{\sigma}.
\end{aligned}
$$

We can identify

$$ dW_{\sigma,u} = \mathbf{u}(x,t)d\mathbf{F}_{\sigma}. \quad (79) $$

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

## Appendix C: Displacement of a Transported Control Surface

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

$$ 
\begin{aligned}
dV(t) &= d \int_{\Omega(t)} 1 \, d\mathbf{x} = \int_{\Omega(t)} \nabla \cdot (\mathbf{u}^* \, d\mathbf{t} + \boldsymbol{\sigma}_t \, d\mathbf{B}_t) \, d\mathbf{x} \\
&= \int_{\delta\Omega(t)} \frac{(\mathbf{u}^* \, d\mathbf{t} + \boldsymbol{\sigma}_t \, d\mathbf{B}_t) \cdot \mathbf{n}}{d\mathbf{x}_{d,t}} \, dS.
\end{aligned}
\quad (80) $$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[20] R. Tailleux. Thermodynamics/dynamics coupling in weakly compressible turbulent stratified fluids. *International Scholarly Research Notices*, 2012, 2012.
References

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_3.png" /></figure>
[21] F. L. Tucciarone, E. Mémin, and L. Li. Primitive equations under location uncertainty: Analytical description and model development. *Stochastic Transport in Upper Ocean Dynamics*, page 287, 2023.

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

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

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

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

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