conjugate, namely

$$ \sigma_t(\mathbf{x})d\mathbf{B}_t = \mathrm{Re}\{ik_{\sigma}^{\perp}\alpha \exp(ik_{\sigma} \cdot \mathbf{x}) d\beta_t\}, \quad \mathbf{a} = |\alpha|^2 k_{\sigma}^{\perp} (k_{\sigma}^{\perp})^T, \quad (16) $$

where $\alpha = \hat{\varphi}(k_\sigma)$ is assumed to be deterministic and real. This monochromatic noise can be directly extended to a multi-scale noise model (15a). Results with a multiscale version of the noise will be shown in the numerical section.

Taking now the expectation $(\mathbb{E})$ of the linearized random system (14), we have

$$ \partial_t \mathbb{E}[\mathbf{u}] + f_0 \mathbb{E}[\mathbf{u}]^{\perp} + g \nabla \mathbb{E}[\eta] - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \mathbb{E}[\mathbf{u}]) = 0, \quad (17a) $$

$$ \partial_t \mathbb{E}[\eta] + H \nabla \cdot \mathbb{E}[\mathbf{u}] - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \mathbb{E}[\eta]) = 0, \quad (17b) $$

where $\partial_t$ denotes the partial time derivative. In order to infer the mean wave solutions, we look for a deterministic ansatz of the form

$$ \mathbb{E}[\tilde{\mathbf{q}}](\mathbf{x}, t) = \hat{\mathbf{q}}_0 \exp (i(\mathbf{k} \cdot \mathbf{x} - \omega t)), \quad (18) $$

where $\mathbf{q} = [\mathbf{u}, \eta]^T = \mathrm{Re}\{\tilde{\mathbf{q}}\}$ is a compact notation for the prognostic variables of the RSW-LU (14), $\hat{\mathbf{q}}_0$ is the initial constant vector (also assumed to be deterministic) and $\omega$ is the time-frequency. We remark that $\mathbb{E}[\tilde{\mathbf{q}}] = \mathbb{E}[\mathbf{q}]$ due to the commutativity of expectation with linear operators.

Injecting next the previous ansatz together with the constant matrix $\mathbf{a}$ (16) into the system (17), we get $L\hat{\mathbf{q}}_0 = 0$ with

$$ L = \left( -i\omega + \frac{1}{2}|\alpha|^2 (\mathbf{k}_\sigma \times \mathbf{k})^2 \right) I_3 + \begin{bmatrix} 0 & -f_0 i g k \\ f_0 & 0 & i g \ell \\ i H k & i H \ell & 0 \end{bmatrix}, \quad (19) $$

where $I_3$ stands for the 3 × 3 identity matrix, $\mathbf{k}_\sigma \times \mathbf{k} = k_\sigma \ell - \ell_\sigma k$ denotes the cross product between the wavenumber vectors $\mathbf{k}_\sigma$ and $\mathbf{k}$ with $|\mathbf{k}_\sigma \times \mathbf{k}| = |\mathbf{k}_\sigma||\mathbf{k}| \sin(\theta)$ and $\theta$ is the angle between them.

As usual, the dispersion relations are then given by the solution of $\det(L) = 0$, namely

$$ \omega = -\frac{i}{2} |\alpha|^2 (\mathbf{k}_\sigma \times \mathbf{k})^2, \quad \omega = -\frac{i}{2} |\alpha|^2 (\mathbf{k}_\sigma \times \mathbf{k})^2 \pm \sqrt{gH|\mathbf{k}|^2 + f_0^2}. \quad (20) $$

We realize immediately that when the noise has the same direction as the initial wave (i.e. $\mathbf{k}_\sigma \times \mathbf{k} = 0$), these two frequencies correspond to the steady and Poincaré (inertia-gravity) waves of the deterministic system [19, 28]. Conversely, when the noise is not aligned to the initial wave, then the term $-\frac{i}{2}|\alpha|^2(\mathbf{k}_\sigma \times \mathbf{k})^2$ leads to
a damping of the ensemble mean wave. We next construct the mean plane wave solutions associated to the frequencies (20).

### 4.1.1 Mean Poincaré Waves

With the value of the last two (opposite) frequencies in (20) and the associated eigenvector of L (19), one obtains the following polarization relations:

$$ \mathbb{E}[\tilde{q}] = \begin{bmatrix} \frac{\omega k + i f_0 \ell}{H |k|^2} \\ \frac{\omega \ell - i f_0 k}{H |k|^2} \\ 1 \end{bmatrix} \hat{\eta}_0 \exp(i(k \cdot x - \omega t)) \exp\left(-\frac{1}{2} |\alpha|^2 (k_\sigma \times k)^2 t\right). \quad (21) $$

Taking the real part, we finally deduce the ensemble-mean wave solution:

$$ \mathbb{E}[\eta] = \hat{\eta}_0 \cos(k \cdot x - \omega t) \exp\left(-\frac{1}{2} |\alpha|^2 (k_\sigma \times k)^2 t\right), \quad (22a) $$

$$ \mathbb{E}[u] = \mathbb{E}[u_\parallel] \frac{k}{|k|} + \mathbb{E}[u_\perp] \frac{k^\perp}{|k|}, \quad (22b) $$

$$ \mathbb{E}[u_\parallel](x, t) = \frac{\hat{\eta}_0 \omega}{H |k|} \cos(k \cdot x - \omega t) \exp\left(-\frac{1}{2} |\alpha|^2 (k_\sigma \times k)^2 t\right), \quad (22c) $$

$$ \mathbb{E}[u_\perp](x, t) = \frac{\hat{\eta}_0 f_0}{H |k|} \sin(k \cdot x - \omega t) \exp\left(-\frac{1}{2} |\alpha|^2 (k_\sigma \times k)^2 t\right), \quad (22d) $$

where the component $\mathbb{E}[u_\parallel]$ is associated to the mean pressure waves that depends on surface elevation mean, whereas the latter component $\mathbb{E}[u_\perp]$ is associated to mean vorticity waves that are initiated by rotation. Note that in the short waves limit with $\|k\|^2 \gg 1/L_d^2$, the mean Poincaré wave corresponds to the inertia-gravity wave of a shallow water system without rotation. The damping term $\exp(-\frac{1}{2}|\alpha|^2(k_\sigma \times k)^2 t)$ depends on the noise's wavelength and variance. This term is zero when the noise and the wave are colinear (i.e. the angle $\theta$ between $k_\sigma$ and $k$ is zero). For high noise magnitude (and $\theta \neq 0$), the damping occurs in a very fast way. In the long wave limit with $\|k\|^2 \ll 1/L_d^2$, the frequency approaches the Coriolis frequency and the damping term is much less important unless the noise amplitude is very high. Nevertheless, the exponential damping in time remains when the noise and the waves are not aligned.

### 4.1.2 Mean Geostrophic Mode

The polarization for the eigenvalue $\omega = -\frac{1}{2}|\alpha|^2(k_\sigma \times k)^2$ reads
$$ \mathbb{E}[\tilde{q}(\mathbf{x}, t)] = \begin{bmatrix} -i \frac{g}{f_0} \ell \\ i \frac{g}{f_0} k \\ 1 \end{bmatrix} \hat{\eta}_0 \exp (i \mathbf{k} \cdot \mathbf{x}) \exp \left( -\frac{1}{2} |\alpha|^2 (\mathbf{k}_\sigma \times \mathbf{k})^2 t \right). \quad (23) $$

The ensemble-mean of the wave solutions are given by

$$ \mathbb{E}[\eta] = f_0 \hat{\eta}_0 \cos (\mathbf{k} \cdot \mathbf{x}) \exp \left( -\frac{1}{2} |\alpha|^2 (\mathbf{k}_\sigma \times \mathbf{k})^2 t \right), \quad (24a) $$

$$ \mathbb{E}[u] = -g \hat{\eta}_0 \mathbf{k}^{\perp} \sin (\mathbf{k} \cdot \mathbf{x}) \exp \left( -\frac{1}{2} |\alpha|^2 (\mathbf{k}_{\sigma} \times \mathbf{k})^2 t \right). \quad (24b) $$

This is a pure vorticity wave. When the noise and the wave are aligned, it corresponds to a steady wave in geostrophic balance, which is also, as we saw, a stationary solution of the nonlinear stochastic system associated to a divergence-free quadratic variation (10). We next look at the path-wise wave solutions under specific noise.

## 4.2 Path-Wise Waves Under Constant Noise

As an initial informative example, we first assume that the noise is constant in space using a zeroth order approximation of the Fourier mode $\exp(ik_\sigma \cdot x)$ in (16). It can be expressed as

$$ \sigma_t d\mathbf{B}_t = \mathrm{Re}\{i\alpha \mathbf{k}_\sigma^\perp d\beta_t\}, \quad \mathbf{a} = \alpha^2 \mathbf{k}_\sigma^\perp (\mathbf{k}_\sigma^\perp)^T. \quad (25) $$

In order to infer wave solutions, we look for stochastic ansatz of the form

$$ \tilde{q}(\mathbf{x}, t) = \hat{q}_0 \exp \left( i (\mathbf{k} \cdot \mathbf{x} - \omega t - \mathrm{Re}\{i \gamma \beta_t\}) \right), \quad (26a) $$

where $\gamma$ is of unit $s^{-1/2}$. This ansatz has been shown to be a solution of a linear stochastic water waves for constant noise in [8]. Applying now the Itô formula [22] for the deterministic function $(t, z) \mapsto \hat{q}_0 \exp (i(c - \omega t - z))$ composed with the random process $\mathrm{Re}\{i\gamma\beta_t\}$, we deduce

$$ d_t \tilde{q} = - \left( \left(i\omega + \frac{1}{2} |\gamma|^2\right) dt + i \mathrm{Re}\{i\gamma d\beta_t\} \right) \tilde{q}, \quad (26b) $$

where the second term on the right-hand-side (RHS) comes from the quadratic variation of the random phase.
Injecting next the stochastic ansatz (26) as well as the noise definition (25) into the linearized system (14), we obtain a system composed of differentiable terms and Brownian (martingale) terms that can be compactly written as

$$L\hat{q}_0 dt = 0, \quad L = \left(-i\omega - \frac{1}{2}|\gamma|^2 + \frac{1}{2}\alpha^2(\mathbf{k}_\sigma \times \mathbf{k})^2\right) I_3 + \begin{bmatrix} 0 & -f_0 & igk \\ f_0 & 0 & ig\ell \\ iHk & iH\ell & 0 \end{bmatrix}, \quad (27a)$$

$$-i\mathrm{Re}\{i\gamma d\beta_t\} + i\mathbf{k} \cdot \mathrm{Re}\{i\alpha \mathbf{k}_\sigma^\perp d\beta_t\} = 0. \quad (27b)$$

The last equation leads to

$$\gamma = \alpha \mathbf{k}_\sigma \times \mathbf{k} \in \mathbb{R}. \quad (28)$$

Substituting it into (27a), we deduce

$$L = \begin{bmatrix} -i\omega - f_0 & igk \\ f_0 & -i\omega & ig\ell \\ iHk & iH\ell & -i\omega \end{bmatrix}. \quad (29a)$$

We remark that in this case the linear operator L for the resolved variables reduces to the same as that of the deterministic system [19, 28]. Solving subsequently det(L) = 0 gives us

$$\omega = 0, \quad \omega = \pm\sqrt{gH|\mathbf{k}|^2 + f_0^2}. \quad (29b)$$

We next construct the stochastic plane wave solutions associated to these frequencies.

### 4.2.1 Stochastic Poincaré Waves

With the value of the last two frequencies in (29b) and the associated eigenvector of L, one can find the following polarization relations:

$$\tilde{q}(\mathbf{x}, t) = \begin{bmatrix} \frac{\omega k + i f_0 \ell}{H |\mathbf{k}|^2} \\ \frac{\omega \ell - i f_0 k}{H |\mathbf{k}|^2} \\ 1 \end{bmatrix} \hat{\eta}_0 \exp \left( i \left( \mathbf{k} \cdot \mathbf{x} - \omega t + \gamma \mathrm{Im}\{\beta_t\} \right) \right). \quad (30)$$

Taking the real part of this ansatz, we deduce the path-wise wave solution:
$$ \eta(\mathbf{x}, t) = \hat{\eta}_0 \cos (\mathbf{k} \cdot \mathbf{x} - \omega t + \gamma \operatorname{Im}\{\beta_t\}), \quad (31a) $$

$$ \begin{aligned} \mathbf{u}(\mathbf{x}, t) = & \frac{\hat{\eta}_0 \omega}{H|\mathbf{k}|} \cos (\mathbf{k} \cdot \mathbf{x} - \omega t + \gamma \operatorname{Im}\{\beta_t\}) \frac{\mathbf{k}}{|\mathbf{k}|} \\ & + \frac{\hat{\eta}_0 f_0}{H|\mathbf{k}|} \sin (\mathbf{k} \cdot \mathbf{x} - \omega t + \gamma \operatorname{Im}\{\beta_t\}) \frac{\mathbf{k}^{\perp}}{|\mathbf{k}|}. \end{aligned} \quad (31b) $$

In this simple case, we can analytically compute the ensemble-mean from these path-wise wave solutions. Note that $X_t := \mathbf{k} \cdot \mathbf{x} - \omega t + \gamma \operatorname{Im}\{\beta_t\} \sim \mathcal{N}(\mathbf{k} \cdot \mathbf{x} - \omega t, \gamma^2 t)$, and the characteristic function of the Gaussian process $X_t$ is given by $\mathbb{E}[\exp(iX_t)] = \exp(i(\mathbf{k} \cdot \mathbf{x} - \omega t) - \frac{1}{2}\gamma^2 t)$. One can then deduce the mean of the random ansatz (26), taking its real part leads finally to

$$ \mathbb{E}[\eta] = \hat{\eta}_0 \cos (\mathbf{k} \cdot \mathbf{x} - \omega t) \exp \left( -\frac{1}{2} \gamma^2 t \right), \quad (32a) $$

$$ \mathbb{E}[\mathbf{u}] = \frac{\hat{\eta}_0}{H|\mathbf{k}|^2} \left( \omega \mathbf{k} \cos (\mathbf{k} \cdot \mathbf{x} - \omega t) + f_0 \mathbf{k}^{\perp} \sin (\mathbf{k} \cdot \mathbf{x} - \omega t) \right) \exp \left( -\frac{1}{2} \gamma^2 t \right). \quad (32b) $$

It can be readily observed from the random dispersion relation (28), that we recover the general mean solution (22) presented in the previous section.

### 4.2.2 Stochastic Geostrophic Mode

The polarization for the eigenvalue $\omega = 0$ reads

$$ \tilde{\mathbf{q}}(\mathbf{x}, t) = \begin{bmatrix} -i \frac{g}{f_0} \ell \\ i \frac{g}{f_0} k \\ 1 \end{bmatrix} \hat{\eta}_0 \exp \left( i \left( \mathbf{k} \cdot \mathbf{x} + \gamma \operatorname{Im}\{\beta_t\} \right) \right). \quad (33) $$

The path-wise wave solutions are then given by

$$ \eta = f_0 \hat{\eta}_0 \cos (\mathbf{k} \cdot \mathbf{x} + \gamma \operatorname{Im}\{\beta_t\}), \quad \mathbf{u} = -g \hat{\eta}_0 \mathbf{k}^{\perp} \sin (\mathbf{k} \cdot \mathbf{x} + \gamma \operatorname{Im}\{\beta_t\}). \quad (34) $$

Similarly, one can recover the general mean solution (24) by taking the expectation of these path-wise solutions. As in the previous case only the waves that are excited by the stochastic forcing remain active on long terms horizon.
## 4.3 Approximation of Path-Wise Waves Under Homogeneous Noise

We now extend the previous solution to statistically homogeneous noise. In the same way as previously we will assume a monochromatic noise as defined in (16), but now slowly varying in space:

$$ \sigma_t dB_t = \text{Re}\{i\alpha k_{\sigma}^{\perp} \exp(i\epsilon k_{\sigma} \cdot x) d\beta_t\}, \quad a = |\alpha|^2 k_{\sigma}^{\perp} (k_{\sigma}^{\perp})^T, \quad (35) $$

where $\epsilon \ll 1$ is a small parameter to ensure the smooth structure of the noise. To infer wave solutions for such homogeneous noise, we now look for the following ansatz:

$$ \tilde{q}(x, t) = \hat{q}_0 \exp \left( i \left( k \cdot x - \omega t - \text{Re}\{i\gamma \exp(i\epsilon k_{\sigma} \cdot x) \beta_t\} \right) \right), \quad (36a) $$

which generalizes our previous ansatz to homogeneous noise. Applying the Itô formula for this ansatz, we have

$$ d_t \tilde{q} = - \left( \left(i\omega + \frac{1}{2}|\gamma|^2\right) dt + i \text{Re}\{i\gamma \exp(i\epsilon k_{\sigma} \cdot x) d\beta_t\} \right) \tilde{q}. \quad (36b) $$

Injecting these solutions ansatz into system (14), we can separate again the Brownian parts and the time-differentiable component. The former reads

$$ -i \text{Re}\{i\gamma \exp(i\epsilon k_{\sigma} \cdot x) d\beta_t\} + i k \cdot \text{Re}\{i\alpha k_{\sigma}^{\perp} \exp(i\epsilon k_{\sigma} \cdot x) d\beta_t\} = 0, \quad (37a) $$

which leads to

$$ \gamma = \alpha k_{\sigma} \times k \in \mathbb{R}. \quad (37b) $$

In a similar way to the previous case, substituting this random dispersion into the diagonal component of the resolved linear operator satisfying $L\hat{q}_0 dt = 0$, we get $\text{diag}(L) = -i\omega I_3$. However, in order to compute the gradient terms of the anti-symmetric part of L, the random phase is linearized as

$$ \begin{aligned} \tilde{q}(x, t) &\approx \hat{q}_0 \exp \left( i \left( k \cdot x - \omega t - \gamma \text{Re}\{i(1 + i\epsilon k_{\sigma} \cdot x)\beta_t\} \right) \right), \\ &= \hat{q}_0 \exp \left( i \left( \underbrace{(k + \epsilon k_{\sigma} \gamma \text{Re}\{\beta_t\})}_{:=\tilde{k}_t} \cdot x - \omega t + \gamma \text{Im}\{\beta_t\} \right) \right). \end{aligned} \quad (38) $$

Hereafter, $\tilde{k}_t = [\tilde{k}, \tilde{\ell}]^T$ is referred to as the effective wavenumber vector ensuing from the space varying random phase. It can be noticed that the real component of the complex Brownian path influences the wave's spatial phase, whereas the wave's temporal phase is randomized by the imaginary component. This latter has already
been considered in the constant noise case. With the previous approximation, the linear operator L can be finally written as

$$L = \begin{bmatrix} -i\omega & -f_0 & ig\tilde{k}_t \\ f_0 & -i\omega & ig\tilde{l}_t \\ iH\tilde{k}_t & iH\tilde{l}_t & -i\omega \end{bmatrix}. \qquad (39a)$$

The two resulting dispersion relations are now given by

$$\omega = 0, \quad \omega = \pm\sqrt{gH|\tilde{k}_t|^2 + f_0^2}. \qquad (39b)$$

The latter random dispersion relation reduces to the previous relation (29b) (associated to a spatially constant noise) when $\epsilon = 0$. We note that the homogeneous random noise leads to wave scattering. Such phenomena corresponds to similar results found in the setting of the Wentzel-Kramers-Brillouin approximation [15, 21] or more recently through Wigner transform [7, 13, 14]. The stochastic framework explored here leads nevertheless to simpler formal developments.

In the same way as previously, we exhibit in the following the two types of stochastic waves associated to this spatially slowly varying homogeneous noise.

### 4.3.1 Stochastic Poincaré Waves

The path-wise wave solution in this case can be written as

$$\eta(\mathbf{x}, t) = \hat{\eta}_0 \cos(\tilde{k}_t \cdot \mathbf{x} - \omega t + \gamma \operatorname{Im}\{\beta_t\}), \qquad (40a)$$

$$\begin{aligned} u(\mathbf{x}, t) = & \frac{\hat{\eta}\omega}{H|\tilde{k}_t|} \cos(\tilde{k}_t \cdot \mathbf{x} - \omega t + \gamma \operatorname{Im}\{\beta_t\}) \frac{\tilde{k}_t}{|\tilde{k}_t|} \\ & + \frac{\hat{\eta}f_0}{H|\tilde{k}_t|} \sin(\tilde{k}_t \cdot \mathbf{x} - \omega t + \gamma \operatorname{Im}\{\beta_t\}) \frac{\tilde{k}_t^{\perp}}{|\tilde{k}_t|}. \end{aligned} \qquad (40b)$$

It can be remarked that the surface elevation phase has two sources of randomness: a temporal one, in the modified frequency, and one in space coming from the space varying noise. For $k = k_σ$ the solutions corresponds again to the classical deterministic waves, while, as shown in Sect. 4.1, when $k \neq k_σ$, the mean of the stochastic wave solutions (22) is damped compared to the deterministic surface elevation ($h_d(\mathbf{x}, t)$).
### 4.3.2 Stochastic Geostrophic Mode

The path-wise wave solutions are given by

$$ \eta = f_0 \hat{\eta}_0 \cos (\tilde{k}_t \cdot \mathbf{x} + \gamma \operatorname{Im}\{\beta_t\}), \quad \mathbf{u} = -g \hat{\eta}_0 \tilde{k}_t^{\perp} \sin (\tilde{k}_t \cdot \mathbf{x} + \gamma \operatorname{Im}\{\beta_t\}). \quad (41) $$

These solutions correspond to dispersive geostrophic modes wave packet. The velocity wave is a pure vorticity wave packet. Its ensemble mean is also damped for $k \neq k_\sigma$ and corresponds to the geostrophic stationary wave for $k = k_\sigma$ (24).

As a final word, on the general linear stochastic shallow water system, it can be noticed that in the long wave limit with $|\tilde{k}|^2 \ll 1/L_d^2$, (when the frequency approaches the Coriolis frequency), the pressure gradient force becomes negligible compared to the other terms in (14a), and for a noise amplitude of order unity the linear shallow water system boils down to a linear stochastic transport equation in a rotating frame:

$$ d_t u + \left( f_0 u^{\perp} - \frac{1}{2} \nabla \cdot (a \nabla u) \right) dt + \sigma_t dB_t \cdot \nabla u = 0. \quad (42) $$

Up to the diffusion and noise term (whose energy balances exactly) this closely corresponds to the so-called near inertial regime in which the fluid is purely transported. In the LU setting, the noise acts always as a random deviation whose energy is exactly compensated by the diffusion term. In the particular case of the shallow water model (linear or nonlinear) the total energy is conserved. For the linear system the total energy of the mean being damped up to a constant, the total energy variance increases up to a constant as a consequence of the total energy conservation.

## 4.4 Numerical Illustrations

We next give simple illustrations of the stochastic wave solutions. Here, rather than evaluating the analytical solutions presented in the previous sections, we propose to discretize numerically the linearized RSW-LU system (14) and perform Monte-Carlo simulations in order to estimate both path-wise and ensemble-mean wave solutions. To that end, the spectral (Fourier) method is adopted for the spatial discretization within a periodic domain, and an exponential integrator [8] combined with the Milstein scheme [9] is used to approximate the mild solution. This semi-discrete problem can be written as

$$ \hat{q}^{(1)} = -ik \cdot (q_t \hat{\sigma}_t \Delta B_t), \quad (43a) $$

$$ \hat{q}^{(2)} = -ik \cdot (q^{(1)} \hat{\sigma}_t \Delta B_t), \quad (43b) $$
$$ \hat{\boldsymbol{q}}_{t+\Delta t} = \exp(A\Delta t) \left( \hat{\boldsymbol{q}}_t + \hat{\boldsymbol{q}}^{(1)} + \frac{1}{2}\hat{\boldsymbol{q}}^{(2)} \right), \quad (43c) $$

$$ A = \begin{bmatrix} 0 & f_0 & -igk \\ -f_0 & 0 & -ig\ell \\ -iHk & -iH\ell & 0 \end{bmatrix}, \quad \boldsymbol{q} = \begin{bmatrix} u \\ v \\ \eta \end{bmatrix}, \quad (43d) $$

where ⋆ denotes the projection coefficient on the discrete Fourier modes, $\boldsymbol{q} = \mathcal{F}^{-1}(\hat{\boldsymbol{q}})$ is the inverse discrete Fourier transform of $\hat{\boldsymbol{q}}$, $\Delta t$ and $\Delta B_t$ stands for the timestep and the Brownian motion's increment respectively. We remark that the classical 2/3 dealiasing rule can be adopted for the practical computations of (43a) and (43b).

A deterministic monochromatic wave corresponding to a single frequency of the Poincaré waves (propagating to the left) is fixed as the initial condition:

$$ \hat{\boldsymbol{q}}_0 = \begin{bmatrix} \frac{\omega k + i f_0 \ell}{H |\boldsymbol{k}|^2} \\ \frac{\omega \ell - i f_0 k}{H |\boldsymbol{k}|^2} \\ 1 \end{bmatrix} \delta(\boldsymbol{k} - \boldsymbol{k}_0), \quad \omega = \sqrt{g H |\boldsymbol{k}|^2 + f_0^2}, \quad (44) $$

where $\delta$ denotes the Dirac function and $\boldsymbol{k}_0$ is the initial wavenumber vector.

For the simulation configuration, we consider a squared shallow basin of length $L = 5120$ km and depth $H = 100$ m at mid-latitude ($f_0 = 10^{-4} \text{s}^{-1}$) and a large-scale wave with $\boldsymbol{k} = [3\Delta k, 0\Delta \ell]^T$, where $\Delta k = \Delta \ell = 2\pi/L$. Each random system either under constant (25) or under homogeneous noises (16) has been simulated over 5 years with 100 realizations. We remark that here we do not use the smoothness approximation (35) in the homogeneous case. The noise's amplitude is fixed as $\alpha = \sqrt{\tau g/f_0}$, where $\tau = \Delta x/\sqrt{gH}$ with $g = 9.81 \text{ m s}^{-2}$, $\Delta x = 40$ km and $\Delta t = 5\tau$.

Figure 1 illustrates the evolution of the path-wise surface elevation $\eta$ and of the ensemble-mean solution $\mathbb{E}[\eta]$ under a constant noise which has a different direction (with $\boldsymbol{k}_\sigma = [4\Delta k, 6\Delta \ell]^T$) than that of the wave ($\boldsymbol{k}_0$). In this case, the path-wise solution preserves the magnitude of the initial monochromatic wave while the ensemble-mean wave is damped along time.

Figure 2 demonstrates the results obtained with the homogeneous noise (with the same $\boldsymbol{k}_\sigma$ as in the previous case). In that case, the path-wise wave is dispersive (scattering effect) whereas the mean wave is dissipative.

The solutions related to the homogeneous noise is rougher than for a constant noise. The damping associated to the mean solution is visually similar in both case.

We then diagnosed an energy decomposition with respect to ensemble of different runs in Fig. 3. This decomposition consists of the (ensemble) mean of energy $\mathbb{E}[\int_D \frac{1}{2}(H|\boldsymbol{u}|^2 + g\eta^2) d\boldsymbol{x}]$, the energy of (ensemble) mean $\int_D \frac{1}{2}(H|\mathbb{E}[\boldsymbol{u}]|^2 + g\mathbb{E}[\eta]^2) d\boldsymbol{x}$ and the energy of “eddy” $\mathbb{E}[\int_D \frac{1}{2}(H|\boldsymbol{u} - \mathbb{E}[\boldsymbol{u}]|^2 + g(\eta - \mathbb{E}[\eta])^2) d\boldsymbol{x}]$. Figure 3 shows that for both noise models, the energy of mean waves is quickly
<figure><img src="image_94.png" /><figcaption>Fig. 1 Illustration for the path-wise surface elevation (top) and its ensemble-mean (bottom) with constant noise (with $k_σ \times k \neq 0$) at different time (by columns)</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>Fig. 2 Illustration for the path-wise surface elevation (top) and its ensemble-mean (bottom) with an homogeneous noise (with $k_σ \times k \neq 0$) at different time (by columns)</figcaption></figure>

dissipates along time while the energy of eddy waves increases with the same rate and continuously backscatters the variance to the ensemble. This mechanism ensures that the mean of the random energy is preserved in time. This corresponds well to the characteristic (2d), (14d) of the proposed stochastic transport model.

Figure 3 illustrates also the dissipation rate of the energy of the mean waves in terms of different scales of the noise (equivalently, different angles of directions
<figure><img src="image_96.png" /><figcaption>**Fig. 3** (Top) Time evolution of energy decomposition (w.r.t. ensemble) with constant (left) and homogeneous (right) noises; (Bottom) Comparison of dissipation of ensemble mean for different noise's scales</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>**Fig. 4** Illustrations under multi-scale noise. From left to right: pathwise wave solution after 5 years simulation, the corresponding ensemble-mean and evolution of energy decomposition</figcaption></figure>

between the noise and the wave). The numerical results confirm our analyses in Sect. 4: the larger the angle $\theta$, or the smaller the noise's scale $k_{\sigma}$, the faster the mean waves are damped in both cases.

Furthermore, propagation of the monochromatic wave by a multi-scale noise model have been numerically tested. In particular, we consider a band of wavenumbers for the noise, $k_{\sigma} = \{k_j\}_{j=m,...,M}$ with $|k_m| < |k_j| < |k_M|$, $\forall m < j < M$, satisfying the $-3$ spectrum power law: $\alpha_j^2|k_j^{\perp}|^2 = \alpha_m^2|k_m^{\perp}|^2 (|k_j|/|k_m|)^{-3}$, i.e. $\alpha_j = \alpha_m (|k_j|/|k_m|)^{-5/2}$. Figure 4 shows the results with $k_m = [5\Delta k, 5\Delta \ell]^T$, $\alpha_m = \alpha$ (as mentioned above) and 10 wavenumbers in total with uniform step $\Delta k$. After 5 years of simulation, the pathwise wave is dispersive while the mean wave is dissipative, and both of them are more irregular than the monochromatic noise solutions (see
Fig. 2). We obtain also the consistent conclusion for the ensemble decomposition of the total energy. Moreover, the conversion from energy of mean to energy of eddy in this case is more faster and efficient than that of the monochromatic noise model (see Fig. 3).

# 5 Shallow Water PV Dynamics and Geostrophic Adjustment

The geostrophic adjustment is the process by which the flow and the pressure field tend to mutually adjust at large scale under the influence of earth's rotation. An obvious manifestation of this adjustment corresponds to isobaric wind field and geostrophic current in the ocean. In the classical case of the deterministic linear shallow water model the geostrophic adjustment can be explained in terms of a variational formulation as the state of minimum energy corresponding to a given constant value of the potential vorticity (PV). Within the LU modelling, the geostrophic adjustment can be explained from the nonlinear system. To that end, let us derive the PV equation of the stochastic system (2). Taking first the curl of the momentum equation (2a) under the space-time invariant Coriolis parameter $f_0$, we have

$$ \mathbb{D}_t (\xi + f_0) = -(\xi + f_0) \nabla \cdot \mathbf{u} \, dt + [\nabla^{\perp} \mathbf{u}, \nabla \sigma_t \, dB_t]_F \\ + \frac{1}{2} \sum_{i,j=1,2} \partial_{x_i} \partial_{x_j} (\nabla a_{ij} \times \mathbf{u}) \, dt, \qquad (45a) $$

where $\xi = \nabla \times \mathbf{u} = \partial_x v - \partial_y u$ denotes the relative vorticity and $[A, B]_F = \text{tr}(A^T B)$ stands for the Frobenius inner product of two matrices $A$ and $B$. Applying next the chain rule [24] for the stochastic transport of height (2b), one obtain

$$ \mathbb{D}_t h^{-1} = h^{-1} \nabla \cdot \mathbf{u} \, dt. \qquad (45b) $$

Let us recall the product rules of two stochastic transport equations in the following. In particular, if two arbitrary tracers $\theta$ and $\zeta$ satisfy

$$ \mathbb{D}_t \theta = \Theta \, dt, \quad \mathbb{D}_t \zeta = Z \, dt + \Sigma dB_t, \qquad (46a) $$

where $\Theta$, $Z$ are time-differentiable forcing terms and $\Sigma$ a martingale forcing component, then according to the Itô-integration-by-part formula [22], we have

$$ \mathbb{D}_t (\theta \zeta) = \theta \mathbb{D}_t \zeta + \zeta \mathbb{D}_t \theta - d \left\langle \int_0^t \sigma_s dB_s \cdot \nabla \theta, \int_0^t \Sigma dB_s \right\rangle_t, \qquad (46b) $$
where the last bracket term denotes the quadratic covariation of two martingales [22]. Applying such rule for $\theta = h^{-1}$ and $\zeta = \xi + f_0$ associated with (45), one deduces the stochastic evolution of the PV, $q = (\xi + f_0)/h$, namely

$$ \mathbb{D}_t q = h^{-1}[\nabla^{\perp} u, \nabla \sigma_t dB_t]_F + \frac{1}{2}h^{-1} \sum_{i,j=1,2} \partial_{x_i} \partial_{x_j} (\nabla a_{ij} \times u) dt \\ - d\langle \int_0^t \sigma_s dB_s \cdot \nabla h^{-1}, \int_0^t [\nabla^{\perp} u, \nabla \sigma_s dB_s]_F \rangle_t. \quad (47a) $$

Opposite to the deterministic shallow water case the PV is not transported by the stochastic flow in general and some source/sink terms appear on the right-hand side of this PV equation. These source/sink terms reflect here the action of the small-scale on the non conservation of PV. In the deterministic context, PV is very sensitive to turbulence and subgrid modelling [3]. The same mechanism is at play here. We can nevertheless explore the condition for which PV remains conserved in the setting of a stochastic modeling of the small scales effect. The first and last terms cancel if the large-scale velocity field and the small-scale random component align with each other up to a uniform vector field. The second term trivially cancels if the random field is homogeneous in space (as in that case $a$ becomes a constant matrix). With these two previous conditions (alignment and homogeneity) PV is path-wise conserved. For homogeneous (incompressible) noise the expectation can be written in flux form and the mean PV is globally conserved. For homogeneous field, the PV equation reduces then to

$$ \mathbb{D}_t q = h^{-1}[\nabla^{\perp} u, \nabla \sigma_t dB_t]_F - d\langle \int_0^t \sigma_s dB_s \cdot \nabla h^{-1}, \int_0^t [\nabla^{\perp} u, \nabla \sigma_s dB_s]_F \rangle_t. \quad (47b) $$

The above equation, combined with the previous results on the stochastic linear dynamics, enables revisiting the mechanism of geostrophic adjustment in the presence of a forcing. The corresponding global pictures is as follows. Let us first consider random fluctuations in the ocean generated by wind forcing at large scales. Then, due to the wavelength mechanism described in the previous section, all the waves that are not aligned with this forcing are smoothed out and eventually annihilated following an exponential decay. The only waves that eventually remain are aligned with the wind forcing, which at large scale corresponds essentially to near inertial waves, and the PV dynamics corresponds then to a pure transport. Thus, the system tends then to relax to a balanced state as in the deterministic case. Conversely, we can devise a similar picture in the atmospheric context, in a configuration where the ocean plays the role of a noise for the atmosphere. The atmospheric waves will evolve toward a near inertial wave field aligned with the noise, following the exact same process. As a result, merging these two perspectives provides an interesting ocean/atmosphere coupling mechanisms of auto adjustment.
At small scales, things are more complicated as the ISD has to be taken into account with an isotropization process that is likely less obvious. As a result the forcing terms in the PV equation remain. Furthermore, in that case the smooth spatial structure assumption of the noise imposes a very low noise amplitude. For the study of strong small-scale forcing this assumption as well as the wave ansatz associated to it have to be revisited. This will be the objective of future works.

# 6 Conclusions

This stochastic extension of the shallow water equations highlights the behavior of the LU setting for large-scale representation of flow dynamics. Opposite to classical eddy viscosity models, which introduces a similar damping term on the waves form, here, the waves that are sustained by the noise term keep their full classical expressions. This provides a simple mechanism for waves selection and an interesting simplified model explaining the emergence of near inertial waves in ocean atmosphere systems. In the LU representation of the linearized shallow water system, the ensemble-mean waves that are not excited by a noise term with the same wavelength vanish exponentially fast whereas the others correspond to the classical deterministic wave solutions. The decay rate depends on the noise variance and on the wavenumber (in a quadratic way). The vanishing is therefore all the more fast for waves with small wavelength. The noise acts hence as a Dirac comb on the mean wave field.

**Acknowledgments** The authors acknowledge the support of the ERC EU project 856408-STUOD. The authors would like to thank Louis Thiry for his helpful comments, remarks and stimulating discussions. We also warmly thank Baylor Fox-Kemper and Oana Lang for fruitful discussions and advises on this work. The code to reproduce the numerical results is available at https://github.com/matlong/sw-wave-lu.

# References

[1] 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*, 50(4):983-1003, 01 Apr. 2020.

[2] J. Berner and Coauthors. Stochastic parameterization: Toward a new view of weather and climate models. *Bull. Amer. Meteor. Soc.*, 98:565-588, 2017.

[3] A. S. Bodner and B. Fox-Kemper. A breakdown in potential vorticity estimation delineates the submesoscale-to-turbulence boundary in large eddy simulations. *Journal of Advances in Modeling Earth Systems*, 12(10):e2020MS002049, 2020. https://doi.org/10.1029/2020MS002049.

[4] 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, 2021.

[5] B. Chapron, P. Dérion, E. Mémin, and V. Resseguier. Large-scale flows under location uncertainty: a consistent stochastic framework. *QJIMS*, 144(710):251-260, 2018.
[6] C. Cotter, D. Crisan, D. Holm, W. Pan, and I. Shevchenko. Numerically modeling stochastic lie transport in fluid dynamics. *SIAM J. on Multiscale Modeling and Simulation*, 17(1):192-232, 2019.

[7] E. Danioux and J. Vanneste. Near-inertial-wave scattering by random flows. *Phys. Rev. Fluids*, 1:033701, Jul 2016.

[8] E. Dinvay and E. Mémin. Hamiltonian formulation of the stochastic surface wave problem. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 478(2265):20220050, 2022.

[9] C. Fiorini, P.-M. Boulvard, L. Li, and E. Mémin. A two-step numerical scheme in time for surface quasi geostrophic equations under location uncertainty. In *B. Chapron, D. Crisan, D. Holm, E. Mémin, and A. Radomska, editors, Stochastic Transport in Upper Ocean Dynamics*, pages 57-67, Cham, 2023. Springer International Publishing.

[10] C. E. Franzke, T. J. O'Kane, J. Berner, P. D. Williams, and V. Lucarini. Stochastic climate theory and modeling. *Wiley Interdisciplinary Reviews: Climate Change*, 6(1):63-78, 2015.

[11] G. Gottwald, D. T. Crommelin, and C. E. Franzke. Stochastic climate theory. In *Nonlinear and Stochastic Climate Dynamics*, pages 209-240. Cambridge University Press, 2017.

[12] D. Holm. Variational principles for stochastic fluid dynamics. *Proc. R. Soc. A*, 471(20140963), 2015.

[13] H. A. Kafiabad, M. A. C. Savva, and J. Vanneste. Diffusion of inertia-gravity waves by geostrophic turbulence. *Journal of Fluid Mechanics*, 869:R7, 2019.

[14] H. A. Kafiabad, J. Vanneste, and W. R. Young. Wave-averaged balance: a simple example. *Journal of Fluid Mechanics*, 911:R1, 2021.

[15] E. Kunze. Near-inertial wave propagation in geostrophic shear. *Journal of Physical Oceanography*, 15:544-565, 1985.

[16] O. Lang, D. Crisan, and E. Mémin. Analytical properties for a stochastic rotating shallow water model under location uncertainty. *Journal of Mathematical Fluid Mechanics*, 25(2):29, 2023.

[17] 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.

[18] L. Li, E. Mémin, and G. Tissot. Stochastic parameterization with dynamic mode decomposition. In *Stochastic Transport in Upper Ocean Dynamics proceedings*, volume 10 of *Mathematics of Planet Earth*, pages 179-193. Springer-Verlag, Mar. 2023.

[19] A. J. Majda. *Introduction to PDEs and Waves for the Atmosphere and Ocean*. Courant Lecture Notes Series. Courant Institute of Mathematical Sciences, 2003.

[20] E. Mémin. Fluid flow dynamics under location uncertainty. *Geophys. & Astro. Fluid Dyn.*, 108(2):119-146, 2014.

[21] C. Mooers. Several effects of a baroclinic current on the cross-stream propagation of inertial-internal waves. *Geophysical Fluid Dynamics*, 6(3):245-275, 1975.

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

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

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

[25] V. Resseguier, W. Pan, and B. Fox-Kemper. Data-driven versus self-similar parameterizations for stochastic advection by Lie transport and location uncertainty. *Nonlinear Processes in Geophysics*, 27(2):209-234, 2020.

[26] G. Tissot, A. Cavalieri, and E. Mémin. Stochastic linear modes in a turbulent channel flow. *J. Fluid Mech.*, 912, 2021.

[27] G. Tissot, A. V. Cavalieri, and E. Mémin. Input-output analysis of the stochastic Navier-Stokes equations: application to turbulent channel flow. *Physical Review Fluids*, 8(3):033904, 2023.

[28] G. Vallis. *Atmospheric and Oceanic Fluid Dynamics*. Cambridge University Press, 2017.
**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_98.png" /></figure>
# Analysis of Sea Surface Temperature Variability Using Machine Learning

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

Check for
updates

Said Ouala, Bertrand Chapron, Fabrice Collard, Lucile Gaultier, and Ronan Fablet

**Abstract** Sea surface temperature (SST) is a critical factor in the global climate system and plays a key role in many marine processes. Understanding the variability of SST is therefore important for a range of applications, including weather and climate prediction, ocean circulation modeling, and marine resource management. In this study, we use machine learning techniques to analyze SST anomaly (SSTA) data from the Mediterranean Sea over a period of 33 years. The objective is to best explain the temporal variability of the SSTA extremes. These extremes are revealed to be well explained through a non-linear interaction between multi-scale processes. The results contribute to better unveil factors influencing SSTA extremes, and the development of more accurate prediction models.

**Keywords** Sea surface temperature · Machine learning · Stochastic models · Extremes

## 1 Introduction

Sea surface temperature (SST) is a critical parameter in the global climate system [1, 2] and plays a vital role in many marine processes, including ocean circulation, evaporation, and the exchange of heat and moisture between the ocean and atmosphere [3, 4, 5].

In recent years, particular attention has been attracted by marine heat waves, when SST largely exceeds the local expected average values [6, 7]. Extreme SST

<footnote>S. Ouala (✉) · R. Fablet
IMT Atlantique, Lab-STICC, Brest, France
e-mail: said.ouala@imt-atlantique.fr</footnote>

<footnote>B. Chapron
Ifremer, LOPS, Plouzané, France</footnote>

<footnote>F. Collard · L. Gaultier
ODL, Locmaria-Plouzané, France</footnote>
can cause coral bleaching [8, 9], with cascading effects on the entire ecosystem. Additionally, localized events affect the amount of atmospheric moisture available, to impact precipitation patterns and the likelihood of drought or flooding in certain regions [10]. Better uncovering factors contributing to these extreme events is therefore of great importance to help predicting and mitigating their impacts.

The SST dynamics compound many processes that interact across a continuum of spatio-temporal scales. A first-order approximation of such a system was initially introduced by [11, 12]. Hasselmann pioneered a two-scale stochastic decomposition, to represent the interactions between slow and fast variables. In this study, we focus on SSTA data collected in the Mediterranean Sea, and examine the potential of machine learning techniques to derive relevant dynamical models. Focus is given on the seasonal modulation of the SSTA and we wish to unveil factors influencing the temporal variability of SSTA extremes. The proposed analysis builds on Hasselmann’s assumption that the variability of the SSTA can be decomposed into slow and fast components. The slow variables mostly follows the seasonal cycle, while the fast variables are linked to rapid processes, e.g. the wind variability. We thus approximate the probability density function of the SSTA data, using a stochastic differential equation in which the drift function represents the seasonal cycle and the diffusion function represents the envelope of the fast SSTA response.

The paper is organized as follows. We start by introducing the general underlying state space model of the SST anomaly. Rather than directly presenting the stochastic model, we first assume that an underlying deterministic ordinary differential equation (ODE) can represent the non-periodic variability of the SSTA. Considering a phase space reconstruction setting, we use the neural embedding of dynamical systems (NbedDyn) framework [13, 14] for this task. We then discuss the limitations of such a representation, and present the stochastic model. We conclude by summarizing our findings and potential future directions.

## 2 Method

Let us assume the following state-space model

$$ \dot{z}_t = f(z_t) \quad (1) $$

$$ x_t = H(z_t) \quad (2) $$

where $t \in [0, +\infty]$ is time. The variables $z_t \in \mathbb{R}^s$ and $x_t \in \mathbb{R}^n$ represent the state variables and the SST anomaly observations respectively. $f$ and $H_t$ are the dynamical and observation operators. The impact of noise on the dynamics and observation models is omitted for simplicity of the presentation.
## 2.1 Deterministic Model Hypothesis

**The NbedDyn Framework** If we assume that $z_t$ is asymptotic to a limit-set $L \subset \mathbb{R}^s$ and that the observations model is not an embedding [15], The NbedDyn model allows one to jointly derive a geometric reconstruction of the unseen phase space from partial observations and a corresponding dynamical model. For any given operator $\mathcal{H}$ of a deterministic dynamical system, Takens theorem [16] guarantees that such an augmented space exists. However, instead of using a delay embedding, NbedDyn defines a $d_E$-dimensional augmented space with states $\mathbf{u}_t \in \mathbb{R}^{d_E}$ as follows:

$$ \mathbf{u}_t^T = [\mathbf{x}_t^T, \mathbf{y}_t^T] \qquad (3) $$

where $\mathbf{y}_t \in \mathbb{R}^{d_E-r}$ are stated as the latent states and $T$ represents the matrix transpose. They account for the unobserved components of the true state $\mathbf{z}_t$.

The augmented state $\mathbf{u}_t$ is assumed to satisfy the following state space model:

$$ \dot{\mathbf{u}}_t = f_{\theta_1}(\mathbf{u}_t) \qquad (4) $$

$$ \mathbf{x}_t = (\mathbf{G}\mathbf{u}_t) \qquad (5) $$

where $\mathbf{G}$ is a projection matrix that satisfies $\mathbf{x}_t = \mathbf{G}\mathbf{u}_t$. The dynamical operator $f_{\theta_1}$ belongs to a given family of neural network operators parameterized by a parameter vector $\theta_1$. In this work, we follow [14] and use a linear quadratic parameterization of $f_{\theta_1}$. This particular parameterization allows us to guarantee boundedness of the ODE (4) using the Schlegel boundedness theorem [17]. A linear quadratic ODE model can be written as follows:

$$ \dot{\mathbf{u}}_t = f_{\theta_1}(\mathbf{u}_t) = \mathbf{c} + \mathbf{L}\mathbf{u}_t + [\mathbf{u}_t^T \mathbf{Q}^{(1)}\mathbf{u}_t, \dots, \mathbf{u}_t^T \mathbf{Q}^{(d_E)}\mathbf{u}_t]^T \qquad (6) $$

where $\mathbf{c} \in \mathbb{R}^{d_E}$, $\mathbf{L} \in \mathbb{R}^{d_E \times d_E}$ and $\mathbf{Q}^{(i)} = [q_{i,j,k}]_{j,k=1}^{d_E}, i = 1, \dots, d_E$. The above approximate model is shifted according to $\bar{\mathbf{u}}_t = \mathbf{u}_t - \mathbf{m}$ with $\mathbf{m} \in \mathbb{R}^{d_E}$. The approximate dynamical equation of the shifted state can be written as:

$$ \dot{\bar{\mathbf{u}}}_t = \mathbf{d} + \mathbf{A}\bar{\mathbf{u}}_t + [\bar{\mathbf{u}}_t^T \mathbf{Q}^{(1)}\bar{\mathbf{u}}_t, \dots, \bar{\mathbf{u}}_t^T \mathbf{Q}^{(d_E)}\bar{\mathbf{u}}_t]^T \qquad (7) $$

with

$$ \mathbf{d} = \mathbf{c} + \mathbf{L}\mathbf{m} + [\mathbf{m}^T \mathbf{Q}^{(1)}\mathbf{m}, \dots, \mathbf{m}^T \mathbf{Q}^{(s)}\mathbf{m}]^T \qquad (8) $$

and
$$ \mathbf{A} = (a_{ij}) = \left( l_{ij} + \sum_{k=1}^{s} (q_{i,j,k} + q_{i,k,j})m_k \right) \qquad (9) $$

Given an observation time series of size $N + 1$ {$\mathbf{x}_{t0}, \dots, \mathbf{x}_{tN}$}, the training setting comes to jointly learning the model parameters $\theta_1 = \{\mathbf{c}, \mathbf{L}, \mathbf{Q}^1, \mathbf{Q}^2, \dots, \mathbf{Q}^d\mathbf{E}, \mathbf{m}\}$ and the latent states $\mathbf{y}_t$ according to the following constrained optimization problem

$$ \begin{aligned} \hat{\theta}_1, \{\hat{\mathbf{y}}_{t_i}\}_{i=0}^{i=N-1} = & \arg \min_{\theta_1, \{\mathbf{y}_{t_i}\}} \sum_{i=1}^{N} \|\mathbf{x}_{t_i} - \mathbf{G}\Phi_{\theta_1, t_i}(\mathbf{u}_{t_{i-1}})\|^2 \\ & + \lambda_1 \|\mathbf{u}_{t_i} - \Phi_{\theta_1, t_i}(\mathbf{u}_{t_{i-1}})\|^2 \\ & + \lambda_2 C_1 \\ & + \lambda_3 C_2 \end{aligned} \qquad (10) $$

with $\Phi_{\theta_1,t}(\mathbf{u}_{t-1}) = \mathbf{u}_{t-1} + \int_{t-1}^{t} f_{\theta_1}(\mathbf{u}_w) dw$ is the flow of the ODE (6) (in our work, this flow is approximated using a Runge Kutta 4 scheme) and $C_1 = \sum_{i,j,k=1}^{s} \|q_{i,j,k} + q_{i,k,j} + b_{j,i,k} + b_{j,k,i} + b_{k,i,j} + b_{k,j,i}\|^2$ and $C_2 = \sum_{i=1}^{s} \text{Max}(\alpha_i, 0)/\text{Max}(\alpha_i + 1, 0)$ where $\alpha_i, i = 1, \dots, d_E$ the eigenvalues of the matrix $A_s = \frac{1}{2}(\mathbf{A} + \mathbf{A}^T)$. The variables $\lambda_{1,2,3}$ are constant weighting parameters. The first constraint $C_1$ steams from the energy-preserving condition of the quadratic non-linearity. It forces the contribution of the quadratic terms of $f_{\theta_1}$ to the fluctuation energy to sum up to zero. The second constraint, $C_2$, ensures that the eigenvalues of $A_s$ are negative. Satisfying these constraints guarantees that the model $f_{\theta_1}$ is bounded through the existence of a monotonically attracting trapping region that includes the limit-set revealed by the minimization of the forecasting loss. Similarly to the Takens delay embedding technique, the sequence:

$$ R_{t_0, t_N} = \{\hat{\mathbf{u}}_{t_i}^T = [\mathbf{x}_{t_i}^T, \hat{\mathbf{y}}_t^T] \text{ with } t_i = t_0, \dots, t_N\} \qquad (11) $$

represents a geometric reconstruction of the phase space. In addition to this reconstruction, the NbedDyn model can be used to forecast new observations by determining an initial condition of the unobserved component $\mathbf{y}_t$ and performing a numerical integration of the ODE model (6). We infer the initial condition using a minimization of an objective function similar to (10), but only with respect to the latent states $\mathbf{y}_t$. This minimization can be seen as a variational data assimilation problem, with partial observations of the state-space variables and known dynamical and observation models [18].

**Related Works** Related state-of-the-art techniques mainly rely on the reconstruction of a phase space using delay embedding [16]. This includes both traditional parametric and non-parametric modeling techniques [19, 20] as well as recurrent neural networks (RNNs). The latter family of methods includes both simple RNN
parameterizations of dynamical systems, as well as latent space inference techniques that are built on an approximation of a posterior distribution that requires the parameterization of a delay embedding [21, 22, 23].

The interest of the NbedDyn framework in contrast to delay embedding based approaches resides in the fact that we do not exploit either a delay embedding or an explicit modeling of the inference model (i.e., the reconstruction of the latent states given the observed time series). As such, our scheme only involves the selection of the class of ODEs of interest. This model reduces the complexity of the overall scheme to the complexity of the ODE representation and guarantees the consistency of the reconstructed latent states w.r.t. the learnt ODE.

## 2.2 Stochastic Model Hypothesis: The Stochastic NbedDyn

When using phase space reconstruction techniques, one should not forget about the assumptions that this theory is built on. For any embedding to work, we are assuming that the dynamical model in (1) exists and can be represented by an ordinary differential equation [15]. For several realistic applications, this ODE may not exist or can have an extremely large dimension. In geoscience, for instance, the dimension of a state space variable can reach $s \approx O(10^9)$. In these situations, reconstructing such an high-dimensional phase space becomes significantly more challenging. In practice, the model returned by any embedding technique can be complemented by an appropriate closure. The form of this closure term can be deterministic using for example the framework of [24] or stochastic through an appropriate calibration of a noise forcing.

When considering SST anomaly data, after calibration of the neural embedding model, an unpredictable, high frequency residual remains. Based on Hasselmann's idea, we assume this residual component represents the effect of fast-scale processes, e.g. passages of atmospheric and oceanic eddies. To first order, it can be represented as a modulated white noise. Indeed, this residual, shown in Fig. 3, exhibits correlations with the slow-scale SST anomaly data.

To model stochastic SST anomalies, the deterministic NbedDyn model described above is first optimized, and further complemented (6) with a stochastic forcing as follows:

$$ \begin{cases} \dot{\mathbf{u}}_t = f_{\theta_1}(\mathbf{u}_t) + g_{\theta_2}(\mathbf{u}_t)\xi_t \\ \mathbf{x}_t = \mathbf{G}\mathbf{u}_t \end{cases} \qquad (12) $$

with $\xi_t$ is a white noise. We derive the parameters of the model (12), as follows. Given an observation time series of size $N + 1$ {$\mathbf{x}_{t_0}, \dots, \mathbf{x}_{t_N}$}, similarly to the deterministic case, we optimize the diffusion parameters $\theta_2$ to minimize the forecast of the observations. In addition to the diffusion parameters, we also reconstruct a noise realization $\xi^{rec}$ that generates the observations process under (12). Overall,
the optimization problem can be written as follows:

$$
\begin{equation}
\begin{aligned}
\hat{\theta}_2, \{\hat{\xi}_{t_i}^{\text{rec}}\}_{i=0}^{i=N-1} &= \arg \min_{\theta_2, \xi^{\text{rec}}} \sum_{t=1}^{T} \|\mathbf{x}_t - \mathbf{G}\mathbf{\Phi}_{\theta,t}(\mathbf{u}_{t-1}, \xi_{t-1}^{\text{rec}})\|^2 \\
&\quad \text{Subject to } \begin{cases} \mathbf{u}_t = \mathbf{\Phi}_{\theta,t}(\mathbf{u}_{t-1}, \xi^{\text{rec}}) \\ \mathbf{G}\mathbf{u}_t = \mathbf{x}_t \\ \mathbf{R}_{\xi^{\text{rec}}\xi^{\text{rec}}}(\tau) = 0 \text{ for all } \tau \neq 0 \end{cases}
\end{aligned}
\tag{13}
\end{equation}
$$

with {$\hat{\xi}_{t_i}^{\text{rec}}$}$_{i=0}^{i=N-1}$ is the noise realization that minimizes the objective function
in (13) and $\mathbf{\Phi}_{\theta,t}$:

$$
\Phi_{\theta,t}(\mathbf{u}_{t-1}, \xi^{\text{rec}}) = \mathbf{u}_{t-1} + \int_{t-1}^{t} f_{\theta}(\mathbf{u}_{w}) dw + \int_{t-1}^{t} g_{\theta}(\mathbf{u}_{w}) \xi_{w}^{\text{rec}} dw
$$

the solution of the stochastic model. This solution is approximated in this work using an Euler-Maruyama scheme, which makes the model converge to an Ito SDE.
In practice, we use the following regularized optimization problem:

$$
\begin{equation}
\begin{aligned}
\hat{\theta}_2, \hat{\xi}^{\text{rec}} ={}& \arg \min_{\theta_2, \xi^{\text{rec}}} \sum_{t=1}^{T} \|\mathbf{x}_t - \mathbf{G}\mathbf{\Phi}_{\theta,t}(\mathbf{u}_{t-1}, \xi^{\text{rec}})\|^2 \\
& + \lambda_4 C_3 \\
& + \lambda_5 C_4 \\
& + \lambda_6 C_5
\end{aligned}
\tag{14}
\end{equation}
$$

with $C_3 = \|\mathbf{R}_{\xi^{\text{rec}}\xi^{\text{rec}}}(\tau)\|^2$, $C_4 = \text{Var}(\Phi_{\theta,t}(\mathbf{u}_{t-1}, \xi^{\text{samp}}))$, $C_5 = \|\Phi_{\theta,t}(\mathbf{u}_{t-1}, \xi^{\text{rec}}) - \mathbf{E}[\Phi_{\theta,t}(\mathbf{u}_{t-1}, \xi^{\text{samp}})]\|^2$ and $\xi^{\text{samp}}$ is a sampled Gaussian white noise. The variables $\lambda_{4,5,6}$ are constant weighting parameters. The first constraint $C_3$ makes the reconstructed noise path white. The second and third constraints, $C_{4,5}$, ensure that the SDE generalizes to sampled white noises. Specifically $C_3$ makes an ensemble of trajectories generated from sampled white noise close to the trajectory generated from the reconstructed noise and $C_2$ reduces the spread of the ensemble around the trajectory simulated from the reconstructed noise.

After optimization, we can couple the optimization problems (10) and (14) and
calibrate jointly all the model parameters θ₁, θ₂, y_t, ξ^rec. This fine tuning step is not
essential but allows both the drift and diffusion parameters of the model to adapt to
each other.
# 3 Numerical Experiments

## 3.1 Data

Sea Surface Temperature Anomalies (SSTA) in the Mediterranean Sea correspond to the Ligurian Sea at 8.6°E, 43.8°N. The anomalies are computed based on a yearly average of the annual 99th percentile of the SST reanalysis [25, 26]. The time series is made up of daily SST anomaly measurements from 1987 to 2019. We use the daily data from 1987 to 2014 as training data. Figure 3e illustrates the time series. These time series include a seasonal cycle and non-periodic high temperature extremes in the summer.

## 3.2 Analysis of the Deterministic Model

In this first experiment, we investigate if the deterministic neural embedding model is able to model the non-periodic variability of the SSTA extremes. For this purpose, we test 3 models with dimensions of the embedding ranging from 1 to 10.

**Analysis of the Embeddings** The choice of the dimension $d_E$ is linked to the number of independent variables that can be used to model the dynamics using, in our context, a bounded autonomous linear quadratic ODE. We start by studying the direct impact of $d_E$ on the performance of the NbedDyn model. Figure 1 shows the impact of $d_E$ on the training error between the observations and the model simulation. Other criteria could be used (please refer to [13, 14] for a more in depth analysis of this parameters on other case studies), but overall, the training error provides a direct measure of the effectiveness of the embedding dimension in the training phase. The first evaluation of the training error reported in Fig. 1

<figure><img src="image_100.png" /><figcaption>**Fig. 1** Mean training error at convergence. We report the mean training error at convergence of the deterministic NbedDyn model for different dimensions $d_E$ of the embedding. This error is averaged over the training time series, and we highlight here both the mean and standard deviations</figcaption></figure>
corresponds to $d_E$ equal to the dimension of the measurements, i.e. $d_E = 1$. In this experiment, no latent states $\mathbf{y}_t$ are used and the embedding $\mathbf{u}_t = \mathbf{x}_t$. In such situations, the ODE model can not perfectly fit the data. Furthermore, at this particular value of $d_E$, the models are more likely to display a bad asymptotic behavior. As the dimension increases, this training error decreases which confirms better modeling abilities using the NbedDyn model. In the following, we study the models with $d_E = 3, 6$ and $10$.

**Asymptotic Properties of the Models** We evaluate the asymptotic behavior of the deterministic models for $d_E = 3, 6$ and $10$. For this purpose, we run the nbedDyn models for a period of 27 years. The resulting simulation is visualized with respect to the reconstructed phase space (11) of the training data in Fig. 2. Overall, the

<figure><img src="image_101.png" /><figcaption>**Fig. 2** Asymptotic solution of the deterministic models. We visualize the simulation of the deterministic models with respect to the reconstructed phase space. The models with $d_E = 3, 6$ and $10$ are given in figure (a), (b) and (c) respectively. For $d_E > 3$, we project the simulation and the reconstructed phase space into $\mathbb{R}^3$</figcaption></figure>
models are only able to reproduce the seasonal cycle of the SST anomaly data. Other experiments (not shown here) suggest that even a farther increase of the dimension of the embedding does not allow the model to capture the non-periodic behavior of the SST anomaly extremes.

**Analysis of the Training Residuals** To further investigate the asymptotic behavior of the deterministic models, we visualize in Fig. 3 the training residual ${x_{t_i} - G\Phi_{\theta,t_i}(u_{t_{i-1}})}$ with $t_i = t_0, \dots, t_N$. When the dimension of the embedding increases above $d_E = 1$, a qualitative and quantitative change in the residual error is present. This is due to the fact that a two-dimensional ODE (in $\mathbb{R}$) is needed to capture the oscillations of the seasonal cycle of the SSTA. However, when the dimension of the embedding increases above 2, no clear qualitative or quantitative change is present. Furthermore, the residual is much more high frequency than the training SSTA data, which suggests that the errors are due to a missing high frequency scales that can not be modeled using the standard deterministic model.

Based on these considerations, and motivated by Hasselmann's works on stochastic climate models with applications on SST anomaly data, we proposed the stochastic NbedDyn model. In this framework, the SST residual of Fig. 3 is modeled as a stochastic forcing.

## 3.3 Analysis of the Stochastic Model

We focus our analysis on the model with $d_E = 6$. We add a stochastic forcing to the neural embedding model (the parameters of the diffusion function $g_{\theta_2}$ are optimized according to Appendix 1). Figure 4 shows the reconstructed phase space under this new model, as well as a model simulation of 27 years. When compared to the simulations of the deterministic model in Fig. 2, the stochastic model is able to cover the whole reconstructed phase space, including the regions with high temperature extremes. This shows that including the high frequency forcing is crucial for the model to capture the non-periodic behavior of the extremes.

These observations are further illustrated in the simulation example given in Fig. 5. The stochastic model is able to produce an ensemble of SST anomaly trajectories that reproduce the non-periodic variability of the extremes. Furthermore, the trajectories generated from a sampled white noise match the one of the reconstructed noise, which validates the proposed training procedure.

We can also discuss the marginal PDF of the stochastic model and compare it to the one computed from the data in Fig. 6. The PDF of the model is computed over a simulation of 109 years. Overall, the model is able to correctly model the high SST anomalies (in the summer), including the non-periodic extremes that form the tail of the distribution. The negative SST anomalies (in the winter) are not approximated as good as in the summer case. This is due to the fact that the model flattens the PDF in the winter by generating trajectories that have more spread (as highlighted in the ensemble prediction experiment in Fig. 5). We did not investigate this problem
<figure><img src="image_102.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_103.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_104.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_106.png" /><figcaption>(e)</figcaption></figure>

**Fig. 3** Training residual and corresponding training data. We visualize the training residual for $d_E = 1$ in (a), $d_E = 3$ in (b), $d_E = 6$ in (c) and $d_E = 10$ in (d). Corresponding training data are given in (e)
<figure><img src="image_107.png" /><figcaption>**Fig. 5** Ensemble simulation of the stochastic model. We visualize an ensemble simulation of the stochastic model, both in the training (a) and test (b) sets. The simulation in the training set is carried to compare trajectories computed from a sampled noise to the one issued from the reconstructed noise $\xi^{rec}$</figcaption></figure>

<figure><img src="image_108.png" /><figcaption>**Fig. 6** PDF of the data and the stochastic model. We compare the marginal PDF of the stochastic model with the one of the data</figcaption></figure>

within the present study. However, we can make the PDF sharper in the winter by forcing the diffusion of the model to be closer to zero during this season.
# 4 Conclusion

In this work, we examined the potential of machine learning techniques to derive relevant dynamical models of sea surface temperature anomaly data in the Mediterranean Sea. We focused on the seasonal modulation of SSTA extremes and used a neural embedding model to reconstruct the phase space of SSTA data. We then added a stochastic forcing term to account for the missing high frequency variability. Our results contribute to the understanding of the factors influencing SSTA extremes and the development of more accurate prediction models. In particular, the analysis highlights the importance of including these fast high-frequency scales in the modeling of SSTA data.

One potential avenue for future work is to investigate the white noise hypothesis in comparison to other types of stochastic models, such as those based on colored noise or fractional Brownian motion. Furthermore, it would be interesting to apply the methodology to other regions and compare the results to evaluate the local impacts of the fast-scales on the slower ones. The ability of using this model as an emulator and studying its predictive skills with respect to standard ocean data assimilation based systems is also a promising perspective.

Finally, and from a methodological point of view, this work highlights the importance of complementing models that are returned by an embedding methodology. Specifically, and as discussed in Sect. 2.2, in complex applications such as the ones in geosciences, the dimension of the underlying state variables is likely to be huge and defining ways of complementing reduced order models through appropriate closure terms is mandatory in order to capture the variability of the data. Analysing the residual of the model fitting procedure is a natural way to define and optimize this closure terms.

# Appendix 1: Training

The trainable parameters of the deterministic NbedDyn models i.e. the linear quadratic ODE and initial conditions of the latent states are initially sampled from a uniform distribution. The training of all models is carried using the Adam optimizer. We use a varying learning rate (from 0.1 to 0.001) in all the experiments to speed up the training. Regarding the weighting parameters {$\lambda_i$}$_{i=1}^{i=3}$, we set $\lambda_1$ = 1 during all the training. The weights responsible for the boundedness constraints were set at higher values in the beginning of the training i.e. $\lambda_2$ = 100 and $\lambda_3$ = 1000 and then reduced to $\lambda_2$ = 1 when $\lambda_3 C_2$ = 0. The training is stopped using cross-validation. Regarding the stochastic forcing, the parameters of the diffusion are initially sampled from a uniform distribution and the noise path $\xi^{rec}$ is initialized from a standard normal distribution. We use a learning rate of 0.001 and the weighting parameters {$\lambda_i$}$_{i=4}^{i=6}$ to one during all the training. We finished the training
with a fine tuning step, in which all the parameters of the model are optimized jointly
with a learning rate of 0.0001.

Appendix 2: Parameterization of the Diffusion Function

The diffusion function $g_{\theta_2} : \mathbb{R}^{d_E} \rightarrow \mathbb{R}^{d_E \times d_N}$ where $d_N$ is the dimension of the noise. In our experiment, we parameterized this function using a fully connected neural networks with 2 hidden layers with a sigmoid activation and 400 neurones per hidden layer. The dimension of the noise $d_N$ is set to the dimension of the state $\mathbf{u}_t$.

References

[1] Grant R Bigg, TD Jickells, PS Liss, and TJ Osborn. The role of the oceans in climate. *International Journal of Climatology: A journal of the Royal Meteorological Society*, 23(10):1127–1159, 2003.

[2] Song Yang, Zhenning Li, Jin-Yi Yu, Xiaoming Hu, Wenjie Dong, and Shan He. El niño–southern oscillation and its impact in the changing climate. *National Science Review*, 5(6):840–857, 2018.

[3] Gerold Siedler, John Gould, and John Church. *Ocean circulation and climate: observing and modelling the global ocean*. Elsevier, 2001.

[4] Alexander Otto, Friederike EL Otto, Olivier Boucher, John Church, Gabi Hegerl, Piers M Forster, Nathan P Gillett, Jonathan Gregory, Gregory C Johnson, Reto Knutti, et al. Energy budget constraints on climate response. *Nature Geoscience*, 6(6):415–416, 2013.

[5] Kevin E Trenberth and John T Fasullo. An apparent hiatus in global warming? *Earth's Future*, 1(1):19–32, 2013.

[6] Eric CJ Oliver, Simon J Wotherspoon, Matthew A Chamberlain, and Neil J Holbrook. Projected Tasman Sea extremes in sea surface temperature through the twenty-first century. *Journal of Climate*, 27(5):1980–1998, 2014.

[7] Deborah Pardo, Stéphanie Jenouvrier, Henri Weimerskirch, and Christophe Barbraud. Effect of extreme sea surface temperature events on the demography of an age-structured albatross population. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 372(1723):20160143, 2017.

[8] Terry P Hughes, James T Kerry, Mariana Álvarez-Noriega, Jorge G Álvarez-Romero, Kristen D Anderson, Andrew H Baird, Russell C Babcock, Maria Beger, David R Bellwood, Ray Berkelmans, et al. Global warming and recurrent mass bleaching of corals. *Nature*, 543(7645):373–377, 2017.

[9] PC Mohanty, A Kushabaha, RS Mahendra, RK Nayak, BK Sahu, E Rao, and T Sinivasa Kumar. Persistence of marine heat waves for coral bleaching and their spectral characteristics around andaman coral reef. *Environmental Monitoring and Assessment*, 193(8):1–9, 2021.

[10] Wenju Cai, Simon Borlace, Matthieu Lengaigne, Peter Van Rensch, Mat Collins, Gabriel Vecchi, Axel Timmermann, Agus Santoso, Michael J McPhaden, Lixin Wu, et al. Increasing frequency of extreme el niño events due to greenhouse warming. *Nature climate change*, 4(2):111–116, 2014.

[11] Klaus Hasselmann. Stochastic climate models part I. Theory. *Tellus*, 28(6):473–485, 1976.
[12] Claude Frankignoul and Klaus Hasselmann. Stochastic climate models, part II application to sea-surface temperature anomalies and thermocline variability. *Tellus*, 29(4):289-305, 1977.

[13] S Ouala, D Nguyen, L Drumetz, B Chapron, A Pascual, F Collard, L Gaultier, and R Fablet. Learning latent dynamics for partially observed chaotic systems. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 30(10):103121, 2020.

[14] Said Ouala, Steven L Brunton, Bertrand Chapron, Ananda Pascual, Fabrice Collard, Lucile Gaultier, and Ronan Fablet. Bounded nonlinear forecasts of partially observed geophysical systems with physics-constrained deep learning. *Physica D: Nonlinear Phenomena*, page 133630, 2023.

[15] Tim Sauer, James A. Yorke, and Martin Casdagli. Embedology. *Journal of Statistical Physics*, 65(3):579-616, Nov 1991.

[16] Floris Takens. Detecting strange attractors in turbulence. In David Rand and Lai-Sang Young, editors, *Dynamical Systems and Turbulence*, Warwick 1980, pages 366-381, Berlin, Heidelberg, 1981. Springer Berlin Heidelberg.

[17] Michael Schlegel and Bernd R. Noack. On long-term boundedness of Galerkin models. *Journal of Fluid Mechanics*, 765:325-352, 2015.

[18] Peter Lynch and Xiang-Yu Huang. *Initialization*, pages 241–260. Springer Berlin Heidelberg, Berlin, Heidelberg, 2010.

[19] Henry D. I. Abarbanel. *Modeling Chaos*, pages 95–114. Springer New York, New York, NY, 1996.

[20] Johan Paduart, Lieve Lauwers, Jan Swevers, Kris Smolders, Johan Schoukens, and Rik Pintelon. Identification of nonlinear systems using Polynomial Nonlinear State Space models. *Automatica*, 46(4):647-656, April 2010.

[21] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep Residual Learning for Image Recognition. arXiv:1512.03385 [cs], December 2015. arXiv: 1512.03385.

[22] Rahul G. Krishnan, Uri Shalit, and David Sontag. Structured Inference Networks for Nonlinear State Space Models. arXiv:1609.09869 [cs, stat], September 2016. arXiv: 1609.09869.

[23] Tian Qi Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. In *Advances in Neural Information Processing Systems*, pages 6571–6583, 2018.

[24] Matthew Levine and Andrew Stuart. A framework for machine learning of model error in dynamical systems. *Communications of the American Mathematical Society*, 2(07):283-344, 2022.

[25] S Simoncelli, C Fratianni, N Pinardi, A Grandi, M Drudi, P Oddo, and S Dobricic. Mediterranean sea physical reanalysis (cmems med-physics)[data set]. *Copernicus Monitoring Environment Marine Service (CMEMS)*, 2019.

[26] Karina von Schuckmann et al. Copernicus marine service ocean state report, issue 3. *Journal of Operational Oceanography*, 12(sup1):S1-S123, 2019.

**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.
# Data Assimilation: A Dynamic Homotopy-Based Coupling Approach

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

Check for
updates

Sebastian Reich

**Abstract** Homotopy approaches to Bayesian inference have found wide-spread use especially if the Kullback-Leibler divergence between the prior and the posterior distribution is large. Here we extend one of these homotopy approaches to include an underlying stochastic diffusion process. The underlying mathematical problem is closely related to the Schrödinger bridge problem for given marginal distributions. We demonstrate that the proposed homotopy approach provides a computationally tractable approximation to the underlying bridge problem. In particular, our implementation builds upon the widely used ensemble Kalman filter methodology and extends it to Schrödinger bridge problems within the context of sequential data assimilation.

## 1 Introduction

Sequential data assimilation interlaces dynamic processes with intermittent partial state observations in order to provide reliable state estimates and their uncertainties. A wide array of numerical methods have been proposed to tackle this problem computationally. Popular methods include sequential Monte Carlo, variational inference, and various ensemble Kalman filter formulations [5, 8]. These methods can encounter difficulties whenever the predictive distribution is incompatible with the incoming data; in other words whenever the distance between the prior, as provided by the underlying stochastic process, and the data informed posterior distribution is large. It has long been realized that this challenge can be partially circumvented by altering the underlying stochastic process through appropriate control terms or modified proposal densities [5, 23, 16, 9]. Recently, the connection between devising such control terms and Schrödinger bridge problems [4] has been made explicit [16]. However, Schrödinger bridge problems are notoriously difficult to

<footnote>S. Reich (✉)
University of Potsdam, Institute of Mathematics, Potsdam, Germany
e-mail: sebastian.reich@uni-postdam.de</footnote>
solve numerically. The key contribution of this paper is to provide a computationally tractable (sub-optimal) solution via a novel extension of established homotopy approaches [7, 14]. Similar to related homotopy approaches for purely Bayesian inference, the solution of certain partial differential equations (PDE) is required in order to find the desired control terms [14, 25]. In line with standard ensemble Kalman filter (EnKF) methodologies we approximate these PDEs via a constant gain approximation [21]. There are also alternative approaches to sequential data assimilation or inference which utilize ideas from optimal transportation; see for example [15, 6, 20, 3].

The paper is structured as follows. The mathematical formulation of the data assimilation problem, as considered in this paper, is laid out in Sect. 2. The standard optimal control and Schrödinger bridge approach to data assimilation is briefly summarized in Sect. 3, and the novel control formulation based on an homotopy formulation is introduced in Sect. 4. A practical implementation based on the EnKF methodology is proposed in Sect. 5. A series of increasingly complex data assimilation problems is considered in Sect. 6 in order to demonstrate the feasibility of the proposed methodologies. The paper concludes with some conclusions in Sect. 7. Detailed mathematical derivations can be found in Appendices 1 and 2, respectively.

## 2 Problem Formulation and Background

We consider drift diffusion processes given by a stochastic differential equations (SDE)

$$dX_t = f(X_t)dt + \sqrt{2\sigma}dW_t, \quad (1)$$

where $X_t : \Omega \to \mathbb{R}^{d_x}$, $f : \mathbb{R}^{d_x} \to \mathbb{R}^{d_x}$, $\sigma \in \mathbb{R}_{\ge 0}$, and $W_t : \Omega \to \mathbb{R}^{d_x}$ denotes $d_x$-dimensional standard Brownian motion [12, 13].

Assuming the law of $X_t$ is absolutely continuous w.r.t. Lebesgue measure with density $\pi_t$, this leads to the Fokker-Planck equation [13]

$$\partial_t \pi_t = -\nabla \cdot (\pi_t (f - \sigma \nabla \log \pi_t)). \quad (2)$$

The SDE (1) can be replaced by the mean field ODE

$$\frac{d}{dt} \tilde{X}_t = f(\tilde{X}_t) - \sigma \nabla \log \tilde{\pi}_t \quad (3)$$

where $\tilde{\pi}_t$ denotes the law of $\tilde{X}_t$. Provided $\tilde{\pi}_0 = \pi_0$, it holds that $\tilde{\pi}_t = \pi_t$ for all $t > 0$. Note that the evolution of the random variable $\tilde{X}_t$ is entirely deterministic subject to random initial conditions $\tilde{X}_0 \sim \pi_0$.

At time $t = T > 0$, we have observations of the system according to
$$y_T = h(x_T^\dagger) + v \quad (4)$$

from which we wish to infer the unknown state $x_T^\dagger \in \mathbb{R}^{d_x}$. Here $h: \mathbb{R}^{d_x} \to \mathbb{R}^{d_y}$ denotes the forward map and $v \sim \mathcal{N}(0, R)$ is $d_y$-dimensional Gaussian noise with covariance matrix $R \in \mathbb{R}^{d_y \times d_y}$.

Let $L: \mathbb{R}^{d_x} \to \mathbb{R}$ denote the corresponding negative log-likelihood function. Since $v$ is Gaussian it is given by

$$L(x) = \frac{1}{2}(h(x) - y_T)^\top R^{-1} (h(x) - y_T) \quad (5)$$

up to an irrelevant constant. The observations are combined with the predictive density $\pi_T$ at time $t = T$ according to Bayes' theorem,

$$\pi_T^a = \frac{e^{-L}\pi_T}{\int e^{-L(x)}\pi_T(x)dx}. \quad (6)$$

The process of transforming the random variable $X_T \sim \pi_T$ into a random variable $X_T^a \sim \pi_T^a$ is called data assimilation in the context of dynamical systems and stochastic processes [10, 17, 8].

Since performing data assimilation can be difficult if the relative Kullback-Leibler divergence

$$\mathrm{KL}(\pi_T | \pi_T^a) = \int_{\mathbb{R}^{d_x}} \pi_T(x) (\log \pi_T(x) - \log \pi_T^a(x)) dx, \quad (7)$$

also called the relative entropy [13], between the prior $\pi_T$ and posterior $\pi_T^a$ is large and/or if the involved distributions are strongly non-Gaussian [19, 1], we propose to construct a new SDE with state process $X_t^h$ such that $X_0^h \sim \pi_0$ and $X_T^h \sim \pi_T^a$. In other words, we are looking for a stochastic process (bridge) with initial density $\pi_0$ and final density $\pi_T^a$. The problem of finding the optimal process (in the sense of minimal Kullback-Leibler divergence) is known as the Schrödinger bridge problem [4].

## 3 Schrödinger Bridge Approach

The Bayesian adjustment (6) at final time $t = T$ leads in fact to an adjustment over the whole solution space of the underlying diffusion process described by (1). Let us denote the so called smoothing distribution by $\pi_t^a$, $t \in [0, T]$ [18, 5]. It is well established that these marginal distributions can be generated from a controlled SDE

$$dX_t^a = f(X_t^a)dt + g_t^a(X_t^a)dt + \sqrt{2\sigma}dW_t \quad (8)$$
for appropriate control $g_t^a : \mathbb{R}^{d_x} \to \mathbb{R}^{d_x}$ such that $X_0^a \sim \pi_0^a$ implies $X_t^a \sim \pi_t^a$ for all $t > 0$. It is also well known that finding a suitable $g_t^a$ can be formulated as an optimal control problem which in turn is closely related to the backward Kolmogorov equation [13, 16]. Formulations related to (8) have also been used in the context of sequential Monte Carlo methods [5].

As proposed in [16], an alternative perspective on sequential data assimilation is provided by Schrödinger bridges. Given two marginal distributions $q_0$ and $q_T$ and stochastic process $X_t$ (referred to as the reference process), a Schrödinger bridge is another stochastic process $\hat{X}_t$ such that $\hat{X}_0 \sim q_0$, $\hat{X}_T \sim q_T$ and the Kullback-Leibler divergence between the processes $\{\hat{X}_t\}_{t \in [0,T]}$ and $\{X_t\}_{t \in [0,T]}$ is minimal. Specialised to our problem and considering a single data assimilation cycle that means the marginals are the initial and posterior densities, i.e. $q_0 = \pi_0$ and $q_T = \pi_T^a$, and the reference process is the solution to (1). The solution to the associated Schrödinger bridge problem is again of the form (8) with modified control term denoted by $g_t^{SB}(x)$.

A Schrödinger bridge is thus the optimal coupling as measured by the Kullback-Leibler divergence to the underlying reference process. Unfortunately Schrödinger bridges lead to boundary value problems in the space of probability measures and the required control term $g_t^{SB}$ seems rather difficult to compute in practice. In addition to the computational complexity of solving nonlinear Schrödinger bridge problems, the target distribution $\pi_T^a$ is implicitly defined in the setting of data assimilation. The next section offers a solution to both of these issues. We point to [23] for a discussion of alternative approaches which introduce appropriate control terms into data assimilation procedures.

# 4 Homotopy Induced Dynamic Coupling

Since Schrödinger bridges are computationally challenging, we ask whether a less optimal but cheaper approach might also be feasible. Indeed, in the context of data assimilation a non-optimal coupling can be found via a homotopy between the initial and target distribution as follows. Let

$$ \pi_t^h(x) = Z_t^{-1} e^{-\frac{t}{T} L(x)} \pi_t(x) \quad (9) $$

denote the homotopy in question, with $Z_t = \int e^{-\frac{t}{T} L(x)} \pi_t(x) dx$ the time dependent normalization constant. It clearly holds that $\pi_0^h = \pi_0$ and $\pi_T^h = \pi_T^a$. Note that the scaling $t \mapsto e^{-\frac{t}{T} L}$ was chosen for its simplicity and follows previous work on Bayesian inference problems [7, 14]. Finding better homotopies or systematic ways of constructing one could be an interesting direction for future research.

We can then reason backwards from the Fokker-Planck equation of $\pi_t^h$ to conclude that if it is the density of a random variable $X_t^h$ then that random variable must satisfy the modified SDE:
$$dX_t^h = f(X_t^h)dt - \frac{\sigma t}{T} \nabla L(X_t^h)dt + g_t(X_t^h)dt + \sqrt{2\sigma} dW_t, \quad (10)$$

where $g_t$ is a solution to the PDE

$$\nabla \cdot (\pi_t^h g_t) = \frac{1}{T} \pi_t^h \left( L + t \nabla L \cdot \left( f - \frac{\sigma t}{T} \nabla L - \sigma \nabla \log \pi_t^h \right) \right) + \pi_t^h \frac{\dot{Z}_t}{Z_t}. \quad (11)$$

The derivations of (11) can be found in Appendix 1.

Note that (10) constitutes a mean field model since $g_t$ depends on the distribution $\pi_t^h$ of $X_t^h$. We also wish to point out that

$$\hat{g}_t^{\text{SB}}(x) := -\frac{\sigma t}{T} \nabla L(x) + g_t(x) \quad (12)$$

provides a solution to the associated coupling problem. The control term $\hat{g}_t^{\text{SB}}$ is however non-optimal in the sense of the Schrödinger bridge problem since it does not minimise the Kullback-Leibler divergence.

Since (11) is linear in $g_t$ we can decompose (11) into a set of simpler equations $\nabla \cdot (\pi_t^h g_t^i) = \pi_t^h (k^i - \mathbb{E}k^i)$ such that the $k^i$ add up to the right hand side of (11). In order to maintain $\int \nabla \cdot (\pi_t^h g_t) = 0$ for the individual $g_t^i$ we can make use of the fact that the terms in (11) are of the form $\pi_t^h (k - \mathbb{E}k)$. Separating the terms, we obtain the following equations, the sum of whose solutions solves (11):

$$\nabla \cdot (\pi_t^h g_t^1) = \pi_t^h \left( \frac{L}{T} - \mathbb{E} \frac{L}{T} \right) \quad (13a)$$

$$\nabla \cdot (\pi_t^h g_t^2) = \pi_t^h \left( \frac{t}{T} \nabla L \cdot f - \mathbb{E} \frac{t}{T} \nabla L \cdot f \right) \quad (13b)$$

$$\nabla \cdot (\pi_t^h g_t^3) = -\pi_t^h \left( \frac{\sigma t}{T} \nabla L \cdot \nabla \log \pi_t^h - \mathbb{E} \frac{\sigma t}{T} \nabla L \cdot \nabla \log \pi_t^h \right) \quad (13c)$$

$$\nabla \cdot (\pi_t^h g_t^4) = -\pi_t^h \left( \frac{\sigma t^2}{T^2} \nabla L \cdot \nabla L - \mathbb{E} \frac{\sigma t^2}{T^2} \nabla L \cdot \nabla L \right). \quad (13d)$$

Note that

$$\pi_t^h \nabla L \cdot \nabla \log \pi_t^h = \nabla (\pi_t^h \nabla L) - \pi_t^h \Delta L, \quad (14)$$

which can be used to avoid the computation of $\nabla \log \pi_t^h$ (with $\Delta = \nabla \cdot \nabla$ the Laplacian operator) in (13c). Thus the controlled SDE (10) can be replaced by

$$dX_t^h = f(X_t^h)dt - \frac{2\sigma t}{T} \nabla L(X_t^h)dt + \hat{g}_t(X_t^h)dt + \sqrt{2\sigma} dW_t, \quad (15)$$
where $\hat{g}_t$ is a solution to

$$ \nabla \cdot (\pi_t^h \hat{g}_t) = \frac{1}{T} \pi_t^h \left( L + t \nabla L \cdot \left( f - \frac{\sigma t}{T} \nabla L \right) + \sigma t \Delta L \right) + \pi_t^h \frac{\dot{Z}_t}{Z_t}. \quad (16) $$

Furthermore, if $\Delta L$ is a constant (as a function of $x$) or small in comparison to the other contributions in (16), then (16) simplifies further. In particular, this is the case if the forward map is linear, that is, $h(x) = Hx$.

Building upon the mean field ODE (3), one obtains the equivalent controlled mean field ODE system

$$ \frac{d}{dt} \tilde{X}_t^h = f(\tilde{X}_t^h) - \frac{\sigma t}{T} \nabla L(\tilde{X}_t^h) + g_t(\tilde{X}_t^h) - \sigma \nabla \log \tilde{\pi}_t^h(\tilde{X}_t^h), \quad (17) $$

with $g_t$ defined as before. This mean field formulation again requires knowledge (or approximation) of $\nabla \log \tilde{\pi}_t^h$. A Gaussian approximation might be sufficient in certain circumstances giving rise to

$$ \nabla \log \tilde{\pi}_t^h(x) \approx -(\tilde{\Sigma}_t^h)^{-1}(x - \tilde{\mu}_t^h), \quad (18) $$

where $\tilde{\mu}_t^h$ denotes the mean of $\tilde{X}_t^h$ and $\tilde{\Sigma}_t^h$ its covariance matrix.

# 5 Numerical Implementation

No analytic solution to (11) is known and we thus have to resort to approximations. We note that a similar PDE arises in the computation of the gain in the feedback particle filter [24] and one could use the diffusion map based approximation [22] for the problem at hand. This method also transforms the PDE into a Poisson equation which it then translated into an equivalent integral equation, the semi-group form of the Poisson equation. As the name suggests the integral equation makes use of the generator of a semi-group which can be approximated by diffusion maps. Here we instead propose to follow the constant gain approximation first introduced in the EnKF methodology [21].

## 5.1 Ensemble Kalman Mean Field Approximation

Let us assume that $\Delta L \approx \text{const}$ in (16). Then we only need to deal with the modified negative log likelihood function

$$ \tilde{L}(x) = \frac{1}{T} \left( L(x) + t \nabla L(x) \cdot \left( f(x) - \frac{\sigma t}{T} \nabla L(x) \right) \right) \quad (19a) $$
$$ \approx \frac{1}{T} \left( L(x) + \frac{t}{\Delta t} \left\{ L(x) - L \left( x - \Delta t f(x) + \Delta t \frac{t \sigma}{T} \nabla L(x) \right) \right\} \right) \quad (19b) $$

with $\Delta t$ being the time-step also used later for time-stepping the evolution equations (10) or (15), respectively. Since $L$ is given by (5), we define the modified forward map

$$ \tilde{h}(x) = h \left( x - \Delta t f(x) + \Delta t \frac{\sigma t}{T} \nabla L(x) \right) \quad (20) $$

and thus

$$ \tilde{L}(x) \approx \frac{t + \Delta t}{2\Delta t T} (h(x) - y_T)^{\top} R^{-1} (h(x) - y_T) - \frac{t}{2\Delta t T} (\tilde{h}(x) - y_T)^{\top} R^{-1} (\tilde{h}(x) - y_T). \quad (21) $$

Following the standard EnKF methodology for quadratic loss functions, this suggests to approximate the drift function $\hat{g}_t$ in (15) as follows:

$$ \hat{g}_t^{\text{KF}}(x) = -\frac{t + \Delta t}{\Delta t T} \Sigma_t^{xh} R^{-1} \left( \frac{1}{2} (h(x) + \pi_t^h[h]) - y_T \right) \quad (22a) $$

$$ + \frac{t}{\Delta t T} \Sigma_t^{x\tilde{h}} R^{-1} \left( \frac{1}{2} (\tilde{h}(x) + \pi_t^h[\tilde{h}]) - y_T \right). \quad (22b) $$

Here we have introduced the notation $\pi_t^h[h]$ to denote the expectation value $\mathbb{E}[h]$ of a function $h(x)$ under the PDF $\pi_t^h$. Furthermore, $\Sigma_t^{xh}$ denotes the correlation matrix between $x$ and $h(x)$ under the PDF $\pi_t^h$ etc. The derivation of (22) can be found in Appendix 2.

## 5.2 Particle Approximation and Time-Stepping

The controlled mean field equations (15) can be implemented numerically by the standard Monte Carlo Ansatz, that is, $M$ particles $X_t^{(i)}$ are propagated according to

$$ dX_t^{(i)} = f(X_t^{(i)})dt - \frac{2\sigma t}{T} \nabla L(X_t^{(i)})dt + \hat{g}_t^{\text{KF}}(X_t^{(i)})dt + \sqrt{2\sigma} dW_t^{(i)} \quad (23) $$

for $i = 1, \dots, M$. The required expectation values in $\hat{g}_t^{\text{KF}}$ are evaluated with respect to the empirical measure

$$ \hat{\pi}_t^h(x) = \frac{1}{2} \sum_{i=1}^{M} \delta(x - X_t^{(i)}). \quad (24) $$
The interacting particle system can be time-stepped using an appropriate adaptation of (61) from Appendix 2. The computation of gradients can be avoided by applying the statistical linearisation (60).

# 6 Examples

We now discuss a sequence of increasingly complex examples. The purpose is both to illuminate certain aspects of the proposed control terms as well as to indicate the computational advantages of the proposed methodology. All examples will be based on linear forward maps $h(x) = Hx$ and, therefore $\Delta L$ is constant and can be ignored.

## 6.1 Pure Diffusion Processes

We set the drift $f$ to zero in (1) and also assuming Gaussian initial conditions. Then the control term (22) gives rise to the mean-field SDE

$$ 
\begin{aligned}
dX_t^h = \sqrt{2\sigma} dW_t - \frac{2\sigma t}{T} H^\top R^{-1} (HX_t^h - y_T) dt & \quad (25a) \\
& - \Sigma_t^h H^\top \left\{ \frac{1}{T} R^{-1} - \frac{2\sigma t^2}{T^2} R^{-1} H H^\top R^{-1} \right\} \left( \frac{1}{2} \left( HX_t^h + H\mu_t^h \right) - y_T \right) dt & \quad (25b)
\end{aligned}
$$

in the limit $\Delta t \to 0$. We note that $X_t^h \sim \pi_t^h$ will remain Gaussian for all times and we denote the mean by $\mu_t^h$ and the covariance matrix by $\Sigma_t^h$. Hence, it holds that $\Sigma_t^{xx} = \Sigma_t^h$ and $\pi_t^h[x] = \mu_t^h$.

Please note that the additional drift term in (25a) is pulling $X_t^h$ towards the observation $y_T$ regardless of the value of $\Sigma_t^h$. It should also be noted that the drift term in (25b) can be both attractive or repulsive with regard to the observation $y_T$ depending on the eigenvalues of

$$ \Omega_t = \frac{1}{T} R^{-1} - \frac{2\sigma t^2}{T^2} R^{-1} H H^\top R^{-1}. \quad (26) $$

The strength of this drift term is moderated by the covariance matrix $\Sigma_t^h$.

We consider a one-dimensional problem with $R = 0.01$, $\sigma = 1$, $H = 1$, $y_T = 1$ and $T = 1$. The initial conditions are Gaussian with mean $\mu_0 = 0$ and variance $\Sigma_0 = 1$. It follows that $\pi_1$ is Gaussian with mean $\mu_1 = 0$ and variance $\Sigma_1 = 2$ and the resulting Gaussian posterior $\pi_1^q$ has mean and variance given by
**Fig. 1** Time evolution of the mean $\mu_t^h$ and the variance $\Sigma_t^h$ under the mean field equations (25). Their values at final agree with the posterior values provided by (27)

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

$$ \mu_1^a = Ky_T \approx 0.9524, \quad \mu_1^a = 2 - 2K \approx 0.0952 \quad (27) $$

with Kalman gain $K = 2/(2 + 0.01) \approx 0.9524$.

In Fig. 1 one can find the time evolution of the mean and the variance under the mean field equations (25). The early impact of the data driven control term on the dynamics is perhaps surprising and quite opposite to the standard sequential approach to data assimilation where one first propagates to final time and only then adjusts according to the available data. It is also worth noticing that the sign of the corresponding $\Omega_t$ changes sign at $t_c = \sqrt{2}/20$ implying that the drift term in (25) has a destabilizing effect on the dynamics for $t > t_c$.

## 6.2 Purely Deterministic Processes

We now set $\sigma = 0$ in (1). We obtain from (22) the mean field ODE system

$$ \frac{d}{dt} X_t^h = f(X_t^h) - \frac{t + \Delta t}{\Delta t T} \Sigma_t^h H^\top R^{-1} \left( \frac{1}{2} H (X_t^h + \mu_t^h) - y_T \right) \quad (28a) $$

$$ + \frac{t}{\Delta t T} \Sigma_t^{x\tilde{h}} R^{-1} \left( \frac{1}{2} (\tilde{h}(X_t^h) + \pi_t^h[\tilde{h}]) - y_T \right) \quad (28b) $$

with

$$ \tilde{h}(x) = Hx - \Delta t H f(x). \quad (29) $$

These equations can be expanded giving rise to

$$ \frac{d}{dt} X_t^h = f(X_t^h) - \frac{1}{T} \left\{ \Sigma_t^h + t \Sigma_t^{xf} \right\} H^\top R^{-1} \left( \frac{1}{2} H (X_t^h + \mu_t^h) - y_T \right) \quad (30a) $$
$$ -\frac{t}{2T} \Sigma_t^h H^\top R^{-1} H \left( f(X_t^h) + \pi_t^h[f] \right) \quad (30b) $$

upon ignoring terms of order $O(\Delta t)$. Unless the drift function $f$ is linear, these mean field equations provide only an approximation to the controlled mean field equations (15).

## 6.3 Linear Gaussian Case

It is instructive to investigate the linear case

$$ f(x) = Fx + b \quad (31) $$

in more detail where again everything remains Gaussian provided $X_0^h$ is Gaussian distributed, that is $\pi_0(x) = N(x; \mu_0, \Sigma_0)$. Under these conditions the densities $\pi_t$ and $\pi_t^h$ will also be Gaussian, we write $\pi_t^h(x) = N(x; \mu_t^h, \Sigma_t^h)$. The associated mean field equations follow from Appendix 2 and are given by

$$ \frac{d}{dt} X_t^h = F X_t^h + b + \sigma (\Sigma_t^h)^{-1} (X_t^h - \mu_t^h) - \frac{2\sigma t}{T} H^\top R^{-1} (H X_t^h - y_T) \quad (32a) $$

$$ -C_t H^\top R^{-1} \left( \frac{1}{2} H \left( X_t^h + \mu_t^h \right) - y_T \right) \quad (32b) $$

$$ -\frac{t}{T} \Sigma_t^h H^\top R^{-1} H \left( \frac{1}{2} F \left( X_t^h + \mu_t^h \right) + b \right). \quad (32c) $$

with

$$ C_t = \frac{1}{T} \Sigma_t^h + \frac{t}{T} \Sigma_t^h F^\top - \frac{2\sigma t^2}{T^2} \Sigma_t^h H^\top R^{-1} H. \quad (33) $$

A qualitative discussion can be performed in the scalar case, that is $d_x = 1$, $H = 1$, $\sigma = 1$, $b = 0$, $T = 1$ and $F = \lambda$. One finds that the control terms involving $F$ stabilize the dynamics whenever $\lambda > 0$. This observation is in line with the fact that the data is crucial only if the dynamics in $X_t$ is unstable, that is, $\lambda > 0$.

We consider a two dimensional diffusion process with state variable $x = (x_1, x_2)^\top$ and linear drift term (31) given by

$$ F = \begin{pmatrix} -2 & 1 \\ 1 & -2 \end{pmatrix}, \quad (34) $$

$b = 0$, and diffusion constant $\sigma = 0.1I$. The forward operator is $H = (1 \ 0)$ and the variance of the noise $R = 0.01$. The initial distribution was $\pi_0 = N((1, 3), 0.02I)$. The observed value at time $T = 1$ is set to $y_T = 2.5$.
<figure><img src="image_111.png" /><figcaption>Fig. 2 Left panel: Time evolution of the mean in $x_1$ and $x_2$, the two associated variances and the covariance between $x_1$ and $x_2$ under the linear diffusion process. Right panel: Time evolution of the same quantities under the controlled diffusion process</figcaption></figure>

The posterior mean takes values $\mu_1^a \approx 2.25$ and $\mu_2^a \approx 1.50$, while the posterior covariance matrix becomes

$$ \Sigma^a \approx \begin{pmatrix} 0.0086 & 0.0039 \\ 0.0039 & 0.0503 \end{pmatrix}. \qquad (35) $$

Numerical results can be found in Fig. 2. The impact of the control term on the linear diffusion process can clearly be seen and is most prominent on the observed $x_1$ component of the process. The final values of the controlled process agree well with their posterior counterparts.

## 6.4 Nonlinear Diffusion Example

We consider a two-dimensional problem and denote the state variable by $x = (x_1, x_2)^T$. The drift term is given by

$$ f(x) = -\nabla V(x), \qquad V(x) = \frac{\lambda_1}{2} (x_2 - 2 + \beta x_1^2)^2 + \frac{\lambda_2}{2} \left( \frac{x_1^4}{2} - x_1^2 \right) \qquad (36) $$

with parameters $\lambda_1 = 2000$, $\lambda_2 = 5$, and $\beta = 1/5$. The diffusion constant is set to $\sigma = 1$. The choice of the potential $V(x)$ has two effects: (1) there is a relative high barrier for particles to pass from positive to negative $x_1$-values and vice versa; (2) the dynamics stay close to the parabola $x_2 = 2 - \beta x_1^2$.

The initial distribution is obtained by sampling $x_1$ from a Gaussian with mean 1.5 and variance 0.0625. The $x_2$ component is obtained from the relation
<figure><img src="image_112.png" /><figcaption>**Fig. 3** Initial (blue) and final particle positions (red) under the given evolution process together with the posterior approximation provided by the EnKF (yellow). The observed value is also displayed</figcaption></figure>

$$x_2 = 2 - \beta x_1^2. \qquad (37)$$

We observe the first component $x_1$ of the state vector at time $T = 1$ with measurement error variance $R = 0.01$. The observed value is set to $y_T = -1.5$. Due to the tiny observation error the posterior is centred sharply about the observed value. Furthermore, recall that the dynamics is essentially slaved to the parabola $x_2 = 2 - \beta x_1^2$ which makes the inference problem strongly nonlinear.

All particle simulations are run with an ensemble size of $M = 1000$. Essentially identical results are obtained for $M = 100$. Smaller ensemble sizes lead to numerical instabilities.

In Fig. 3, one can find the particle distribution at time $t = 1$ which constitutes the prior distribution for the associated Bayesian inference problem. It is obvious that a particle filter would fail to recover the posterior distribution which is sharply centered about the observed value. We found that increasing the ensemble size to $M = 10,000$ allows a particle filter to recover the posterior distribution; but the effective ensemble size still drops dramatically. The approximation provided by the EnKF is also displayed. The EnKF fails to recover the posterior due to its inherent linear regression ansatz which is inappropriate for this strongly nonlinear inference problem even in the limit ensemble size $M \to \infty$.

In Fig. 4, the results from the controlled mean field formulation are displayed. It can be concluded that the posterior distribution is well approximated despite the constant gain approximation made in order to formulate the control term $\hat{g}_t^{\text{KF}}$ in (22).

## 6.5 Lorenz-63 Example

All examples so far have considered a single data assimilation cycle only. We now perform a proper sequential data assimilation experiment for the standard Lorenz-63 model [11]
<figure><img src="image_113.png" /><figcaption>Fig. 4 Left panel: Initial and final particle positions under the controlled evolution process. Right panel: Particle positions at intermediate times $t_k \in [0, 1]$</figcaption></figure>

$$ \frac{d}{dt} X_t = f(X_t), \qquad (38) $$

where $X_t : \Omega \to \mathbb{R}^3$ and

$$ f(x, y, z) = \begin{pmatrix} a(y - x) \\ x(b - z) - y \\ xy - cz \end{pmatrix} \qquad (39) $$

with parameters $a = 10, b = 28$ and $c = 8/3$.

In order to obtain a reference solution $X_t^\dagger$ for $t \ge 0$, the ODE (38) is solved numerically with step-size $\Delta t = 0.005$ and initial condition

$$ X_0^\dagger = \begin{pmatrix} -0.587276 \\ -0.563678 \\ 16.8708 \end{pmatrix}. \qquad (40) $$

Scalar-valued observations are generated every $\Delta t_{\text{obs}} > 0$ units of time using the forward model

$$ y_n \Delta t_{\text{obs}} = H X_n \Delta t_{\text{obs}} + v_n, \quad n = 1, \dots, N, \qquad (41) $$

with measurement errors $v_n \sim N(0, 1)$ and forward map $H = (100) \in \mathbb{R}^{1 \times 3}$. We use $\Delta t_{\text{obs}} \in \{0.05, 0.1, 0.12\}$ in our experiments and perform $N = 20,000$ assimilation cycles.

The initial ensemble $\{X_0^{(i)}\}_{i=1}^M$ is drawn from the Gaussian distribution with mean $X_0^\dagger$ and covariance matrix $0.01I$. We employ multiplicative ensemble inflation which amounts to replacing the Lorenz-63 dynamics by
$$ \frac{d}{dt} X_t^{(i)} = f(X_t^{(i)}) + \sigma_k (X_t^{(i)} - \hat{\mu}_t), \quad i = 1, \dots, M, \quad (42) $$

with inflation factors

$$ \sigma_k = 0.025k, \quad k = 0, \dots, 9. \quad (43) $$

Here $\hat{\mu}_t$ denotes the empirical mean of the ensemble $\{X_t^{(i)}\}_{i=1}^M$. These equations are combined with the augmented evolution equations (30) and solved numerically with step-size $\Delta t = 0.005$ and ensemble sizes $M \in \{5, 10, 15\}$.

We report the resulting root mean square errors

$$ \text{RMSE} = \sqrt{\frac{1}{3N} \sum_{n=1}^{N} \| \hat{\mu}_{n \Delta t_{\text{obs}}} - X_{n \Delta t_{\text{obs}}}^{\dagger} \|^2}, \quad (44) $$

which are computed for each ensemble size $M$, observation interval $\Delta t_{\text{obs}}$ and inflation factor $\sigma_k$. The results are displayed in Table 1 where the smallest RMSE over the range of inflation factors $\{\sigma_k\}_{k=0}^9$ is stated for each $M$ and $\Delta t_{\text{obs}}$. We also state the corresponding RMSEs from a standard ensemble square root filter implementation [2, 8]. We find that the proposed homotopy approach outperforms the ensemble square root filter in terms of RMSE in all settings considered. The improvements increase for increasing observation intervals $\Delta t_{\text{obs}}$. The homotopy approach also appears less sensitive to the ensemble size $M$.

We close this example by pointing out that less of an improvement could be expected for a fully observed Lorenz-63 system. The proposed homotopy approach seems particularly effective in guiding the unobserved solution components to regions of high posterior probability. See also the example from Sect. 6.4.

<table><caption>**Table 1** RMSE for both a standard ensemble square root filter (ESRF) implementation and our homotopy approach in terms of ensemble sizes $M \in \{5, 10, 15\}$ and observation intervals $\Delta t_{\text{obs}} \in \{0.05, 0.1, 0.12\}$. The homotopy based data assimilation method leads to significantly reduced RMSEs in all settings</caption><thead><tr><th>M/Δt<sub>obs</sub></th><th></th><th>5</th><th>10</th><th>15</th></tr></thead><tbody><tr><td>0.05</td><td>ESRF</td><td>0.57</td><td>0.56</td><td>0.57</td></tr><tr><td>0.05</td><td>Homotopy</td><td>0.55</td><td>0.55</td><td>0.55</td></tr><tr><td>0.10</td><td>ESRF</td><td>0.85</td><td>0.82</td><td>0.82</td></tr><tr><td>0.10</td><td>Homotopy</td><td>0.77</td><td>0.76</td><td>0.77</td></tr><tr><td>0.12</td><td>ESRF</td><td>0.96</td><td>0.95</td><td>0.94</td></tr><tr><td>0.12</td><td>Homotopy</td><td>0.86</td><td>0.86</td><td>0.86</td></tr></tbody></table>
# 7 Conclusions

Devising alternative proposal densities has a long history in the context of sequential data assimilation and filtering. Here we have explored a computationally tractable approach which combines the concept of Schrödinger bridges with a rather straight-forward homotopy approach. A further key ingredient is the approximate solution of the arising PDEs in terms of a constant gain approximation, which is also widely used within the EnKF community. Numerical examples indicate that the approach is viable and can overcome limitations of both standard sequential Monte Carlo as well as standard EnKF methods. This has been demonstrated for single assimilation steps as well as long-time data assimilation using the chaotic Lorenz-63 model with only the first component observed infrequently. It remains to be seen how the proposed methods behave for high dimensional stochastic processes.

**Acknowledgments** This research has been funded by the Deutsche Forschungsgemeinschaft (DFG)- Project-ID 318763901 - SFB1294. We thank Nikolai Zaki for earlier work on the topic of this paper.

# Appendix 1: Derivation of Control Term Equation

Given an evolution equation

$$dX_t = f(X_t)dt + \sqrt{2\sigma}dW_t \quad (45)$$

we obtain the Fokker-Planck equation

$$\partial_t \pi_t = -\nabla \cdot (\pi_t (f - \sigma \nabla \log \pi_t)). \quad (46)$$

Now we modify (45) by an additional drift term, i.e.

$$dX_t^h = f(X_t^h)dt + \tilde{g}_t(X_t^h)dt + \sqrt{2\sigma}dW_t \quad (47)$$

with $\tilde{g}_t: \mathbb{R}^{d_x} \rightarrow \mathbb{R}^{d_x}$. In that case, we would get a Fokker-Planck equation for the new equation:

$$\partial_t \pi_t^h = -\nabla \cdot (\pi_t^h (f + \tilde{g} - \sigma \nabla \log \pi_t^h)) = -\nabla \cdot (\pi_t^h (f - \sigma \nabla \log \pi_t^h)) - \nabla \cdot (\pi_t^h \tilde{g}_t). \quad (48)$$

We can find $\tilde{g}_t$ in terms of known quantities as follows: we begin by taking the derivative of $\pi_t^h$ with respect to time:

$$\partial_t \pi_t^h = -\pi_t^h \left( \frac{\dot{Z}_t}{Z_t} + \frac{L}{T} \right) + \frac{1}{Z_t} e^{-\frac{t}{T}L} \partial_t \pi_t. \quad (49)$$
Next we substitute (46) for $\partial_t \pi_t$ and use $\pi_t = Z_t e^{\frac{t}{T} L} \pi_t^h$:

$$ \partial_t \pi_t^h = - \pi_t^h \left( \frac{\dot{Z}_t}{Z_t} + \frac{L}{T} \right) - Z_t^{-1} e^{-\frac{t}{T} L} \nabla \cdot \left( \pi_t (f - \sigma \nabla \log \pi_t) \right) \quad (50a) $$

$$ = - \pi_t^h \left( \frac{\dot{Z}_t}{Z_t} + \frac{L}{T} \right) - Z_t^{-1} e^{-\frac{t}{T} L} \nabla \cdot \left( Z_t e^{\frac{t}{T} L} \pi_t^h \left( f - \sigma \nabla \log Z_t e^{\frac{t}{T} L} \pi_t^h \right) \right) \quad (50b) $$

$$ = - \pi_t^h \left( \frac{\dot{Z}_t}{Z_t} + \frac{L}{T} \right) - e^{-\frac{t}{T} L} \nabla \cdot \left( e^{\frac{t}{T} L} \pi_t^h \left( f - \sigma \frac{t}{T} \nabla L - \sigma \nabla \log \pi_t^h \right) \right) \quad (50c) $$

$$ = - \pi_t^h \left( \frac{\dot{Z}_t}{Z_t} + \frac{L}{T} \right) - \nabla \cdot \left( \pi_t^h \left( f - \sigma \frac{t}{T} \nabla L - \sigma \nabla \log \pi_t^h \right) \right) \quad (50d) $$

$$ - \frac{t}{T} \pi_t^h \nabla L \cdot \left( f - \sigma \frac{t}{T} \nabla L - \sigma \nabla \log \pi_t^h \right). \quad (50e) $$

Comparing with (48) it follows that we require

$$ \nabla \cdot (\pi_t^h \tilde{g}_t) = \frac{1}{T} \pi_t^h \left( L + t \nabla L \cdot \left( f - \sigma \frac{t}{T} \nabla L - \sigma \nabla \log \pi_t^h \right) \right) + \pi_t^h \frac{\dot{Z}_t}{Z_t} - \nabla \cdot \left( \pi_t^h \sigma \frac{t}{T} \nabla L \right). \quad (51) $$

For the $\dot{Z}_t$ term we have

$$ \frac{\dot{Z}_t}{Z_t} = \frac{1}{Z_t} \int \partial_t e^{-\frac{t}{T} L(x)} \pi_t(x) dx \quad (52a) $$

$$ = \frac{1}{Z_t} \int -\frac{L}{T} e^{-\frac{t}{T} L(x)} \pi_t(x) + e^{-\frac{t}{T} L(x)} \partial_t \pi_t(x) dx \quad (52b) $$

$$ = -\frac{1}{T} \mathbb{E}L - \frac{1}{Z_t} \int e^{-\frac{t}{T} L(x)} \nabla \cdot \left( \pi_t (f - \sigma \nabla \log \pi_t) \right) dx \quad (52c) $$

$$ = -\frac{1}{T} \mathbb{E}L - \frac{1}{Z_t} \int \frac{t}{T} e^{-\frac{t}{T} L(x)} \pi_t \nabla L \cdot \left( f - \sigma \nabla \log \pi_t \right) dx \quad (52d) $$

$$ = -\frac{1}{T} \mathbb{E}L - \frac{t}{T} \mathbb{E} \nabla L \cdot \left( f - \sigma \nabla \log \pi_t \right) \quad (52e) $$

$$ = -\frac{1}{T} \mathbb{E}L - \frac{t}{T} \mathbb{E} \nabla L \cdot \left( f - \sigma \nabla \log \pi_t^h - \sigma \frac{t}{T} \nabla L \right), \quad (52f) $$

where the third equality follows from integration by parts and the expected value is with taken with respect to $\pi_t^h$. We finally note that
$$ \tilde{g}_t(x) = g_t(x) - \sigma \frac{t}{T} \nabla L(x). \qquad (53) $$

# Appendix 2: Ensemble Kalman Filter Approximations

We provide details on the derivation of the EnKF-like approximation (22) to the controlled mean field equation (10) and the various simplifications that arise from assuming a linear forward map.

We first recall that a continuous time formulation of the EnKF for a generic likelihood function

$$ L(x) = \frac{1}{2}(\hat{h}(x) - y_T)^{\top}\hat{R}^{-1}(\hat{h}(x) - y_T) \qquad (54) $$

is provided by

$$ \frac{d}{dt}X_t = -\Sigma_t^{x\hat{h}}\hat{R}^{-1}\left(\frac{1}{2}(\hat{h}(X_t) + \pi_t[h]) - y_T\right). \qquad (55) $$

Here $\pi_t$ denotes the law of $X_t$, $\pi_t[h]$ the expectation value of a function $g$ under $\pi_t$, and

$$ \Sigma_t^{x\hat{h}} = \pi_t\left[(x - \pi_t[x])(\hat{h} - \pi_t[h])^\top\right] \qquad (56) $$

is the covariance matrix between the state $x$ and the forward map $\hat{h}$.

Formal application of this approach to the two contributions to the likelihood function (21) leads to (22). More precisely, the first term leads to $\hat{h} = h$ and

$$ \hat{R} = \frac{t + \Delta t}{2\Delta t T} R \qquad (57) $$

while the second term results in $\hat{h} = \tilde{h}$ and

$$ \hat{R} = -\frac{t}{\Delta t T} R. \qquad (58) $$

The EnKF makes use of statistical linearization

$$ \Sigma_t^{xx} \pi_t[\nabla h] = \Sigma_t^{xh} \qquad (59) $$

which holds provided $\pi_t$ is Gaussian or if $h$ is linear; a result known as Stein's identity. The identity can also be used to approximate derivatives in a (weakly) non-Gaussian setting giving rise to
$$ \nabla h(x) \approx (\Sigma_t^{xx})^{-1} \Sigma_t^{xh}. \quad (60) $$

We also recall the robust time-stepping method

$$ X_{t_{n+1}} - X_{t_n} = -\Delta t \Sigma_{t_n}^{x\hat{h}} (\Delta t \Sigma_{t_n}^{\hat{h}} + R)^{-1} \left( \frac{1}{2} (\hat{h}(X_{t_n}) + \pi_{t_n}[h]) - y_T \right), \quad (61) $$

which again can be adjusted appropriately to (22).

We now assume a linear forward map, that is $h(x) = Hx$, and discuss the simplifications that result in the computation of (22). Note that

$$ \tilde{h}(x) = Hx - \Delta t H f(x) + \Delta t \frac{\sigma t}{T} H H^{\top} R^{-1} H x. \quad (62) $$

Hence the covariance matrix $\Sigma_t^{x\hat{h}}$ can be reformulated to

$$ \Sigma_t^{x\hat{h}} = \Sigma_t^{xx} H^{\top} - \Delta t \Sigma_t^{xf} H^{\top} + \Delta t \frac{\sigma t}{T} \Sigma_t^{xx} H^{\top} R^{-1} H H^{\top} \quad (63) $$

and (22) simplifies to

$$ \hat{g}_t^{\text{KF}}(x) = -\frac{1}{T} \Sigma_t^{xx} H^\top R^{-1} \left( \frac{1}{2} (Hx + H\mu_t^h) - y_T \right) \quad (64\text{a}) $$

$$ -\frac{t}{T} \Sigma_t^{xf} H^\top R^{-1} \left( \frac{1}{2} (Hx + H\mu_t^h) - y_T \right) \quad (64\text{b}) $$

$$ + \frac{\sigma t^2}{T^2} \Sigma_t^{xx} H^\top R^{-1} H H^\top R^{-1} \left( \frac{1}{2} (Hx + H\mu_t^h) - y_T \right) \quad (64\text{c}) $$

$$ -\frac{t}{2T} \Sigma_t^{xx} H^\top R^{-1} H \left( f(x) + \pi_t^h[f] \right) \quad (64\text{d}) $$

$$ + \frac{\sigma t^2}{T^2} \Sigma_t^{xx} H^\top R^{-1} H H^\top R^{-1} \left( \frac{1}{2} (Hx + H\mu_t^h) - y_T \right) + \mathcal{O}(\Delta t) \quad (64\text{e}) $$

$$ = -C_t H^\top R^{-1} \left( \frac{1}{2} (Hx + H\mu_t^h) - y_T \right) \quad (64\text{f}) $$

$$ -\frac{t}{2T} \Sigma_t^{xx} H^\top R^{-1} H \left( f(x) + \pi_t^h[f] \right) + \mathcal{O}(\Delta t) \quad (64\text{g}) $$

with

$$ C_t = \frac{1}{T} \Sigma_t^{xx} + \frac{t}{T} \Sigma_t^{xf} - \frac{2\sigma t^2}{T^2} \Sigma_t^{xx} H^\top R^{-1} H. \quad (65) $$
Upon dropping terms of order $O(\Delta t)$ and using $\Sigma_t^{xx} = \Sigma_t^h$, we obtain (25) for $f=0$ and (30) for $\sigma=0$ as special cases. The mean field equations (32) also follow easily from $f(x) = Fx + b$.

## References

[1] S. Agapiou, O. Papaspiliopoulos, D. Sanz-Alonso, and A. Stuart. Importance sampling: Intrinsic dimension and computational cost. *Statistical Science*, pages 405–431, 2017.

[2] M. Asch, M. Bocquet, and M. Nodet. *Data assimilation: Methods, algorithms, and applications*. SIAM, 2016.

[3] E. Bernton, J. Heng, A. Doucet, and P. E. Jacob. Schrödinger bridge samplers. Technical report, arXiv:1912.13170, 2019.

[4] Y. Chen, T. T. Georgiou, and M. Pavon. Stochastic control liaisons: Richard Sinkhorn meets Gaspard Monge on a Schrödinger bridge. *SIAM Review*, 63:249–313, 2021. doi: https://doi.org/10.1137/20M1339982.

[5] N. Chopin and O. Papaspiliopoulos. *An Introduction to Sequential Monte Carlo*. Springer Nature Switzerland AG, Cham, Switzerland, 2020.

[6] A. Corenflos, J. Thornton, G. Deligiannidis, and A. Doucet. Differentiable particle filtering via entropy-regularized optimal transport. In *International Conference on Machine Learning*, pages 2100–2111. PMLR, 2021.

[7] F. Daum, J. Huang, and A. Noushin. Exact particle flow for nonlinear filters. In I. Kadar, editor, *Signal Processing, Sensor Fusion, and Target Recognition XIX*, volume 7697, pages 92–110. International Society for Optics and Photonics, SPIE, 2010. doi: https://doi.org/10.1117/12.839590.

[8] G. Evensen, F. Vossepoel, and P. van Leeuwen. *Data Assimilation Fundamentals: A unified Formulation of the State and Parameter Estimation Problem*. Springer Nature Switzerland AG, Cham, Switzerland, 2022.

[9] J. Heng, A. N. Bishop, G. Deligiannidis, and A. Doucet. Controlled sequential Monte Carlo. *The Annals of Statistics*, 48:2904–2929, 2020. doi: https://doi.org/10.1214/19-AOS1914.

[10] K. Law, A. Stuart, and K. Zygalakis. Data assimilation. Cham, Switzerland: Springer, 214, 2015.

[11] E. Lorenz. Deterministic non-periodic flows. *J. Atmos. Sci.*, 20:130–141, 1963.

[12] B. Oksendal. *Stochastic differential equations: an introduction with applications*. Springer Science & Business Media, 2013.

[13] G. A. Pavliotis. *Stochastic Processes and Applications: Diffusion Processes, the Fokker-Planck and Langevin Equations*, volume 60 of *Texts in Applied Mathematics*. Springer New York, New York, NY, 2014. doi: https://doi.org/10.1007/978-1-4939-1323-7.

[14] S. Reich. A dynamical systems framework for intermittent data assimilation. *BIT Numerical Mathematics*, 51(1):235–249, 2011.

[15] S. Reich. A nonparametric ensemble transform method for Bayesian inference. *SIAM J. Sci. Comput.*, 35:A2013–A2024, 2013. doi: https://doi.org/10.1137/130907367.

[16] S. Reich. Data assimilation: The Schrödinger perspective. *Acta Numerica*, 28:635–711, 2019.

[17] S. Reich and C. Cotter. *Probabilistic forecasting and Bayesian data assimilation*. Cambridge University Press, 2015.

[18] S. Särkkä. *Bayesian Filtering and Smoothing*. Cambridge University Press, Cambridge, 2013.

[19] C. Snyder, T. Bengtsson, P. Bickel, and J. Anderson. Obstacles to high-dimensional particle filtering. *Monthly Weather Review*, 136(12):4629–4640, 2008.

[20] A. Spantini, R. Baptista, and Y. Marzouk. Coupling techniques for nonlinear ensemble filtering. *arXiv preprint arXiv:1907.00389*, 2019.
[21] A. Taghvaei, J. de Wiljes, P. Mehta, and S. Reich. Kalman filter and its modern extensions for the continuous-time nonlinear filtering problem. *ASME. J. Dyn. Sys., Meas., Control.*, 140: 030904-030904-11, 2017. doi: https://doi.org/10.1115/1.4037780.

[22] A. Taghvaei, P. Mehta, and S. Meyn. Diffusion map-based algorithm for gain function approximation in the feedback particle filter. *SIAM/ASA J. Uncertain. Quantif.*, 8(3):1090-1117, 2020. doi: https://doi.org/10.1137/19M124513X.

[23] P. van Leeuwen, H. Künsch, L. Nerger, R. Potthast, and S. Reich. Particle filters for high-dimensional geoscience applications: A review. *Q. J. Royal Meteorol. Soc.*, 145:2335-2365, 2019. doi: https://doi.org/10.1002/qj.3551.

[24] T. Yang, P. G. Mehta, and S. P. Meyn. Feedback particle filter. *IEEE Trans. Automat. Control*, 58(10):2465-2480, 2013. ISSN 0018-9286. doi: https://doi.org/10.1109/TAC.2013.2258825.

[25] T. Yang, H. A. P. Blom, and P. G. Mehta. The continuous-discrete time feedback particle filter. In *American Control Conference*, pages 648-653. IEEE, 2014. doi: https://doi.org/10.1109/ACC.2014.6859259.

**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_114.png" /></figure>
<figure><img src="image_115.png" /></figure>

# 1 Introduction

For ensemble-based data assimilation purposes, there is a definite need for relevant ensemble sampling tools. Indeed, the quality and spreading of these ensembles have deep implications in the quality of the data assimilation (Dufée et al 2022), and—until recently—those so-called covariance inflation tools have mostly relied on unsuitable linear Gaussian frameworks (Tandeo et al 2020, Resseguier et al 2020a). A promising alternative is the generation of ensembles through a stochastic remapping of the physical space.

Consider a random mapping $T$, acting at every infinitesimal time step, such that $T_t(x) - x$ is interpreted as a “location perturbation” expressed by

$$T_t(x) = x + a(t, x)\Delta t + e_i(t, x)\Delta\eta_i(t), \quad (1.1)$$

where $a(t, x)$, $e_i(t, x) \in \mathbb{R}^n$. In Eq. (1.1), $a(t, x)$ controls deterministic location shifts, and $\Delta\eta_i(t) \sim \mathcal{N}(0, \Delta t)$ random ones. At every time step, this random mapping $T$ shall induce a perturbation to any tensor field $\theta(t)$ (Zhen et al 2022). For instance, one can perturb a differential form $\theta(t)$ applying $\theta(t) \to T_t^*\theta(t)$ with $T_t^*$ the associated pull-back operator.

<footnote>V. Resseguier (✉)
INRAE, OPAALE, Rennes, France
LAB, SCALIAN DS, Rennes, France
e-mail: valentin.resseguier@scalian.com</footnote>

<footnote>Y. Zhen
Department of Oceanography, Hohai University, Nanjing, China</footnote>

<footnote>B. Chapron
Laboratoire d'Océanographie Physique et Spatiale (LOPS), Ifremer, Plouzané, France</footnote>
A rigorous mathematical definition and calculation of $T_t$ and $T_t^*$ can be obtained in terms of stochastic flows of diffeomorphisms and its Lie derivatives (e.g., Bethencourt De Leon 2021). Yet, to rapidly assess $T_t^*\theta$, a Taylor expansion and Itô's lemma can be used. Given coordinates $(x^1, \dots, x^n)$, when $\theta$ is a differential $k$-form, it can be written as

$$ \theta = \sum_{i_1 < \dots < i_k} f^{i_1 \dots i_k} dx^{i_1} \wedge \dots \wedge dx^{i_k}, \qquad (1.2) $$

with $f$ a semimartingale smooth in space. Then

$$ T_t^*\theta = \sum_{i_1 < \dots < i_k} f^{i_1 \dots i_k} (T_t(x)) T_t^*(dx^{i_1} \wedge \dots \wedge dx^{i_k}), \qquad (1.3) $$

leading to a compact expression

$$ T_t^*\theta = \theta + \mathcal{M}(t, \theta)\Delta t + \mathcal{N}_i(t, \theta)\Delta \eta_i(t), \qquad (1.4) $$

with some differential $k$-forms $\mathcal{M}(t, \theta)$ and $\mathcal{N}_i(t, \theta)$. Appendix 4 provides definitions of $\mathcal{M}$ and $\mathcal{N}$ (see Zhen et al 2022, Appendix B for a full proof).

Hereafter, we present and discuss the potential of this random mapping scheme to possibly prescribe $\theta$, and the parameters $a$ and $e_i$ to ensure that certain quantities, i.e. mass, vorticity, helicity, energy, are conserved.

Several examples of $T_t^*\theta$ can indeed be considered. For instance, when $\theta = f$ is a function (differential 0-form),

$$ (T_t^*\theta) = f + \underbrace{(a^j \partial_{x^j} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f)}_{=\mathcal{M}} \Delta t + \underbrace{e_i^p \partial_{x^p} f}_{=\mathcal{N}_i} \Delta \eta_i. \qquad (1.5) $$

And when $\theta = f dx^1 \wedge \dots \wedge dx^n$ (differential $n$-form), it then follows

$$ T_t^*\theta = \left\{ f + \left( (\partial_{x^p} a^p + \frac{1}{2} J_i) f + (a^p + e_i^p \partial_{x^q} e_i^q) \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f \right) \Delta t \right. \\ \left. + (\partial_{x^p} e_i^p f + e_i^p \partial_{x^p} f) \Delta \eta_i \right\} dx^1 \wedge \dots \wedge dx^n. \qquad (1.6) $$

Finally, when $\theta = f^j dx^j = \sum_{j=1}^n f^j dx^j$ is differential 1-form, we have

$$ T_t^*\theta = \left\{ f^j + (a^p \partial_{x^p} f^j + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f^j + \partial_{x^j} a^p f^p + \partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p) \Delta t \right. \\ \left. + (e_i^p \partial_{x^p} f^j + \partial_{x^j} e_i^p f^p) \Delta \eta_i \right\} dx^j. \qquad (1.7) $$
## 2 Induced Stochastic PDE

From the expression of $T_t^*\theta$, a SPDE is derived from an original PDE, when $\theta$ is a differential form. Suppose $S^d$ is the full state variable of the deterministic dynamical system:

$$ \frac{\partial S^d}{\partial t} = g(S^d). \qquad (2.1) $$

Let $f^d$ be a component or a collection of components of $S^d$. We then associate $f^d$ to a differential form $\theta^d$, i.e. there is an invertible map $\mathcal{F}$ that maps the space of $f^d$ to the space of $\theta^d$, such that $\mathcal{F}(f^d) = \theta^d$. Typically, if $f^d$ is a tracer, it is often associated to the 0-form $\theta^d = f^d$. If $f^d$ is the density $\rho^d$, we might associate the $n$-form $\theta^d = \rho^d dx^{i_1} \wedge \cdots \wedge dx^{i_n}$. More generally, $\theta^d$, and thus $\mathcal{F}$, can be prescribed to ensure that certain quantities—such as mass, energy, circulation—are conserved (Zhen et al 2022, section 3.3). Consider the propagation equation for $f^d$

$$ \delta f^d = g^f(S^d) \delta t. \qquad (2.2) $$

It implies a propagation equation for $\theta$:

$$ \delta\theta^d = g^\theta(S^d)\delta t. \qquad (2.3) $$

We will now stochastically perturb the above deterministic dynamics. Let us denote $S$, $f$ and $\theta$ the semimartingale solutions of this randomized dynamics. The proposed discrete-time perturbation at each time step consists of the following two steps:

$$ \tilde{\theta}(t + \Delta t) = \theta(t) + g^\theta(S(t))\Delta t, \qquad (2.4) $$

$$ \theta(t + \Delta t) = T_t^* \tilde{\theta}(t + \Delta t), \qquad (2.5) $$

with $T_t^* \tilde{\theta}(t + \Delta t) = \tilde{\theta}(t + \Delta t) + \mathcal{M}(t, \tilde{\theta}(t + \Delta t)) \Delta t + \mathcal{N}_i(t, \tilde{\theta}(t + \Delta t)) \Delta \eta_i(t) + o(\Delta t)$ for the associated differential forms $\mathcal{M}(t, \tilde{\theta})$ and $\mathcal{N}_i(t, \tilde{\theta})$.

The deterministic PDE (2.4) and $\|\tilde{\theta}(t + \Delta t) - \theta(t)\|$ scales in $O(\Delta t)$. There is no noise term to induce a scaling in $O(\sqrt{\Delta t})$. Therefore, it can be assumed that there exists $C > 0$ so that $\|\mathcal{M}(t, \tilde{\theta}(t + \Delta t)) - \mathcal{M}(t, \theta(t))\| < C\Delta t$ and $\|\mathcal{N}_i(t, \tilde{\theta}(t + \Delta t)) - \mathcal{N}_i(t, \theta(t))\| < C\Delta t$, for $\Delta t$ small enough. Accordingly,

$$ \begin{aligned} T_t^* \tilde{\theta}(t + \Delta t) = & \tilde{\theta}(t + \Delta t) + (\mathcal{M}(t, \theta(t)) + \mathcal{O}(\Delta t)) \Delta t \\ & + (\mathcal{N}_i(t, \theta(t)) + \mathcal{O}(\Delta t)) \Delta \eta_i(t) + o(\Delta t), \\ = & \tilde{\theta}(t + \Delta t) + \mathcal{M}(t, \theta(t)) \Delta t + \mathcal{N}_i(t, \theta(t)) \Delta \eta_i(t) + o(\Delta t). \end{aligned} \qquad (2.6) $$
Therefore,

$$ \theta(t + \Delta t) = \theta(t) + g^{\theta}(S(t))\Delta t + \mathcal{M}(t, \theta(t))\Delta t + \mathcal{N}_i(t, \theta(t))\Delta \eta_i + o(\Delta t). \quad (2.7) $$

It suggests the following stochastic propagation equation for $\theta$:

$$ \ddot{\theta} = g^{\theta}(S)\dot{\theta}t + \mathcal{M}(t, \theta)\dot{\theta}t + \mathcal{N}_i(t, \theta)\dot{\eta}_i. \quad (2.8) $$

Since there is a 1-1 correspondence between $\theta$ and $f$, Eq. (2.3) also suggests a stochastic propagation equation for $f$, which can be written as

$$ \ddot{f} = g^f(S)\dot{f}t + \mathcal{M}^f(f)\dot{f}t + \mathcal{N}_i^f(f)\dot{\eta}_i. \quad (2.9) $$

We denote the additional terms in Eq. (2.9) by

$$ \ddot{\mathcal{N}}_s f := \mathcal{M}^f(f)\dot{f}t + \mathcal{N}_i^f(f)\dot{\eta}_i. \quad (2.10) $$

Then Eq. (2.9) can be written as:

$$ \ddot{f} = g^f(S)\dot{f}t + \ddot{\mathcal{N}}_s f. \quad (2.11) $$

# 3 Comparison with Other Perturbation Schemes

Obtained above, $\ddot{\mathcal{N}}_s f$ is completely determined by $T_t^*\theta$, but is not directly related to the original dynamics Eq. (2.2). Once the expression of $T$ in Eq. (1.1) and the choice of the differential form $\theta^d$ are determined, the perturbation term $\ddot{\mathcal{N}}_s f$ is prescribed. However, the choice of $\theta^d$ is up to the user, and may then be related to the original dynamics.

In the following, we thus demonstrate that both the stochastic advection by Lie transport (SALT) equation (Holm 2015) and the location uncertainty (LU) equation (Mémin 2014, Resseguier et al 2017; 2020b) can be properly recovered using the proposed perturbation scheme.

## 3.1 Comparison with the LU Equations

The Reynolds transport theorem is central to the LU setting. The Reynolds transport theorem expresses an integral conservation equation for the transport of any conserved quantity within a fluid, connected to its corresponding differential equation. A link between the proposed perturbation approach and the LU formulation can
be anticipated to be related to differential $n$-forms. But first, we consider a key ingredient of LU: the stochastic material derivative of functions (differential 0-forms).

### 3.1.1 0-Forms in the LU Framework

Dropping the forcing terms, LU equation for compressible and incompressible flow reads (Resseguier et al 2017)

$$ \partial_t f + \boldsymbol{w}^* \cdot \nabla f = \nabla \cdot (\frac{1}{2} \boldsymbol{a} \nabla f) - \boldsymbol{\sigma} \dot{\boldsymbol{B}} \cdot \nabla f, \quad (3.1) $$

$$ \boldsymbol{w}^* = \boldsymbol{w} - \frac{1}{2}(\nabla \cdot \boldsymbol{a})^{\top} + \boldsymbol{\sigma}(\nabla \cdot \boldsymbol{\sigma})^{\top}, \quad (3.2) $$

where $\boldsymbol{a} = \boldsymbol{\sigma}_{\bullet k} \boldsymbol{\sigma}_{\bullet k}^T$ and $f$ can be any quantity that is assumed to be transported by the flow, i.e. $Df/Dt = 0$ where $D/Dt$ is the Itô material derivative. For instance, $f$ could be the velocity (dropping forces in the SPDE), the temperature, or the buoyancy.

Separating the terms of the SPDE related to the deterministic dynamics from the term associated to the stochastic scheme, it comes

$$ \ddot{\boldsymbol{a}}^{\text{LU}} f = g^f(S) \dot{\boldsymbol{a}}t + \ddot{\boldsymbol{a}}_s^{\text{LU}} f, \quad (3.3) $$

where

$$ g^f(S) = -\boldsymbol{w} \cdot \nabla f, \quad (3.4) $$

$$ \ddot{\boldsymbol{a}}_s^{\text{LU}} f = -(\boldsymbol{w}^* - \boldsymbol{w}) \cdot \nabla f \dot{\boldsymbol{a}}t - \boldsymbol{\sigma} \ddot{\boldsymbol{a}} \boldsymbol{B} \cdot \nabla f + \nabla \cdot (\frac{1}{2} \boldsymbol{a} \nabla f) \dot{\boldsymbol{a}}t. \quad (3.5) $$

Besides, from our proposed scheme applied to a 0-form $\theta = f$ (Eq. (1.5)), we obtain:

$$ \ddot{\boldsymbol{a}}_s f = \left( a^p \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f \right) \dot{\boldsymbol{a}}t + e_i^p \partial_{x^p} f \ddot{\boldsymbol{a}} \eta_i. \quad (3.6) $$

To physically interpret this equation, we rewrite:

$$ \frac{\partial_s f}{\partial t} = -V^p \partial_{x^p} f + \partial_{x^p} \left( \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f \right), \quad (3.7) $$

where

$$ V^p = -a^p + \frac{1}{2} \partial_{x^q} (e_i^p e_i^q) - e_i^p \frac{\partial \eta_i}{\partial t}. \quad (3.8) $$

Terms of advection and diffusion are recognized. The matrix $\frac{1}{2}e_i e_i^T$ is symmetric non-negative and represents a diffusion matrix. The $p$-th component of the advect-
ing velocity $V^p$ is composed of the drift $-a^p$, a correction $\frac{1}{2}\partial_{x^q}(e_i^p e_i^q)$, and a stochastic advecting velocity $-e_i^p \frac{\partial \eta_i}{\partial t}$.

Direct calculation yields that Eq. (3.5) coincides with Eq. (3.7) when $\mathbf{a} = \boldsymbol{\sigma}_{\bullet k} \boldsymbol{\sigma}_{\bullet k}^T = e_i e_i^T$ and $\boldsymbol{\sigma} \dot{\mathbf{B}} = -e_i d\eta_i$ and

$$T_t(x) = x + e_i^q \partial_{x^q} e_i \Delta t + e_i \Delta \eta_i = x - \mathbf{w}_S^c \Delta t + (-\mathbf{w}_S^c \Delta t - \boldsymbol{\sigma} \Delta \mathbf{B}), \quad (3.9)$$

where

$$\mathbf{w}_S^c = -\frac{1}{2}(\boldsymbol{\sigma}_{\bullet k} \cdot \nabla)\boldsymbol{\sigma}_{\bullet k} = -\frac{1}{2}(\nabla \cdot \mathbf{a})^\top + \frac{1}{2}\boldsymbol{\sigma}(\nabla \cdot \boldsymbol{\sigma})^\top. \quad (3.10)$$

The LU equation can thus be derived by choosing $\theta = f$ and $T_t$ by Eq. (3.9). Note, the term $(-\mathbf{w}_S^c \Delta t - \boldsymbol{\sigma} \Delta \mathbf{B}) = (\frac{1}{2}e_i^q \partial_{x^q} e_i \Delta t + e_i \Delta \eta_i)$ is the Itô noise plus its Itô-to-Stratonovich correction. Hence, it corresponds to the Stratonovich noise $e_i \circ d\eta_i$ of the flow associated to $T_t$. The additional drift $-\mathbf{w}_S^c \Delta t$ is different in nature. It is related to the advecting correction $\mathbf{w}_S^c \cdot \nabla f$ in the LU setting. Indeed, in the LU framework, the Itô drift, $\mathbf{w}$, is seen as the resolved large-scale velocity. That is why, in this framework, the deterministic dynamics (3.4) involve the Itô drift, $\mathbf{w}$. This is also the reason why, under the LU derivation, the advected velocity is assumed to be given by the Itô drift, $\mathbf{w}$. It differs from the Stratonovich drift $\mathbf{w}_S = \mathbf{w} + \mathbf{w}_S^c$, used as advected velocity in SALT approach or in Mikulevicius and Rozovskii (2004) (where the Stratonovich drift is denoted $u$). Interested readers are referred to (Resseguier et al 2020b, Appendix A) for a discussion on these assumptions and for the complete table of SALT-LU notations correspondences. Note however that in all these approaches, the advecting velocity is always the Stratonovich drift. This can be seen e.g., in the Stratonovich form of LU equations, derived in (Resseguier 2017, Appendix 10.1) and (Resseguier et al 2020a, 6.1.3):

$$\partial_t f + \mathbf{w}_S \cdot \nabla f = -(\boldsymbol{\sigma} \circ \dot{\mathbf{B}}) \cdot \nabla f, \quad (3.11)$$

where $\boldsymbol{\sigma} \circ \dot{\mathbf{B}}$ is the Stratonovich noise of the SPDE. Since the advecting velocity $\mathbf{w}_S$ and the resolved velocity $\mathbf{w}$ differ by a drift $\mathbf{w}_S^c$, the term $\mathbf{w}_S^c \cdot \nabla f$ is interpreted as an advection correction, being part of the stochastic scheme (3.5). Accordingly, the remapping $T_t$ involves an additional drift $-\mathbf{w}_S^c \Delta t$.

To also understand (3.9), the inverse flow can be considered:

$$T_t^{-1}(x) = x - e_i \Delta \eta_i = x + \boldsymbol{\sigma} \Delta \mathbf{B}. \quad (3.12)$$

If $T_t$ represents a necessary perturbation to match, at each time step, a true solution, $T_t^{-1}$ measures the difference, at each time step, between this true solution and a model forecast. Therefore, the LU equation can be derived using the proposed perturbation scheme, choosing $\theta = f$ and assuming that a true solution differs from a model forecast by a displacement prescribed by Eq. (3.12).
### 3.1.2 n-Forms in the LU Framework

The LU physical justification relies on a stochastic interpretation of fundamental conservation laws, typically conservation of extensive properties (i.e. integrals of functions over a spatial volume) like momentum, mass, matter and energy (Resseguier et al 2017). These extensive properties can be expressed by integrals of differential n-forms. For instance, the mass and the momentum are integrals of the differential n-forms $\rho dx^1 \wedge \cdots \wedge dx^n$ and $\rho w dx^1 \wedge \cdots \wedge dx^n$, respectively. In the LU framework, a stochastic version of the Reynolds transport theorem (Resseguier et al 2017, Eq. (28)) is used to deal with these differential n-forms $\theta = f dx^1 \wedge \cdots \wedge dx^n$. Assuming an integral conservation $\frac{d}{dt} \int_{V(t)} f = 0$ on a spatial domain $V(t)$ transported by the flow, it leads to the following SPDE:

$$ \frac{Df}{Dt} + \nabla \cdot (\mathbf{w}^* + \boldsymbol{\sigma} \dot{\mathbf{B}}) f = \frac{d}{dt} \left\langle \int_0^t D_t f, \int_0^t \nabla \cdot \boldsymbol{\sigma} \dot{\mathbf{B}} \right\rangle = (\nabla \cdot \boldsymbol{\sigma}_{\bullet i})(\nabla \cdot \boldsymbol{\sigma}_{\bullet i})^T f, \quad (3.13) $$

where $D/Dt$ denotes the Itô material derivative. Forcing terms are dropped for the sake of readability. This SPDE can be rewritten using the expression of that material derivative (Eq. (9) and (10) of Resseguier et al (2017)):

$$ \partial_t f + \nabla \cdot (\mathbf{w}_S f) = \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla f) + \frac{1}{2} \nabla \cdot (\boldsymbol{\sigma}_{\bullet i} (\nabla \cdot \boldsymbol{\sigma}_{\bullet i})^T f) - \nabla \cdot (\boldsymbol{\sigma} \dot{\mathbf{B}} f). \quad (3.14) $$

The original deterministic equation and stochastic perturbation correspond to

$$ g^f(S) = -\nabla \cdot (\mathbf{w} f), \quad (3.15) $$

$$ \begin{aligned} \mathbf{d}_s^{\text{LU}} f = & (-\nabla \cdot (\mathbf{w}_S^c f) + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla f) + \frac{1}{2} \nabla \cdot (\boldsymbol{\sigma}_{\bullet i} (\nabla \cdot \boldsymbol{\sigma}_{\bullet i})^T f)) \mathrm{d}t - \nabla \cdot (\boldsymbol{\sigma} \mathbf{d} \mathbf{B} f), \\ = & -\nabla \cdot ((-\frac{1}{2} \nabla \cdot \mathbf{a})^T \mathrm{d}t + \boldsymbol{\sigma} \mathbf{d} \mathbf{B}) f) + \nabla \cdot (\frac{1}{2} \mathbf{a} \nabla f) \mathrm{d}t. \end{aligned} \quad (3.16) $$

We can now compare these LU equations to our new stochastic scheme applied to n-form $\theta = f dx^1 \wedge \cdots \wedge dx^n$ (Eq. (1.6)). This implies that

$$ \begin{aligned} \mathbf{d}_s f = & ((\partial_{x^p} a^p + \frac{1}{2} J_i) f + (a^p + e_i^p \partial_{x^q} e_i^q) \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f) \mathrm{d}t \\ & + (\partial_{x^p} e_i^p f + e_i^p \partial_{x^p} f) \mathrm{d}\eta_i, \end{aligned} \quad (3.18) $$

where $J_i = \partial_{x^p} e_i^p \partial_{x^q} e_i^q - \partial_{x^p} e_i^q \partial_{x^q} e_i^p$. Rewritten, it leads to:

$$ \frac{\mathrm{d}_s f}{\mathrm{d}t} = -\partial_{x^p} (\tilde{V}^p f) + \partial_{x^p} \left( \left( \frac{1}{2} e_i^p e_i^q \right) \partial_{x^q} f \right), \quad (3.19) $$
where

$$ \tilde{V}^p = V^p - (e_i^p \partial_{xq} e_i^q) = -a^p + \frac{1}{2}(\partial_{xq} e_i^p e_i^q - e_i^p \partial_{xq} e_i^q) - e_i^p \frac{\partial \eta_i}{\partial t}. \quad (3.20) $$

Again an advection-diffusion equation is recognized, but of different nature. Indeed, as expected for an n-form, the PDE is similar to a density conservation equation. Moreover, the advecting drift is slightly different to take into account the cross-correlations between $f(T_t(x))$ and $T_t^*(dx^1 \wedge \cdots \wedge dx^n)$.

Identifying $\mathbf{a} = \boldsymbol{\sigma}_{\bullet k} \boldsymbol{\sigma}_{\bullet k}^T = e_i e_i^T$ and $\boldsymbol{\sigma} \dot{\mathbf{B}} = -e_i d\eta_i$,

$$ \tilde{V} = -a^p + \frac{1}{2}(\partial_{xq} e_i^p e_i^q - e_i^p \partial_{xq} e_i^q) - e_i^p \frac{\partial \eta_i}{\partial t} = -\left(\frac{1}{2} \nabla \cdot \mathbf{a}\right)^T + \boldsymbol{\sigma} \dot{\mathbf{B}}, \quad (3.21) $$

i.e.

$$ a^p = \frac{1}{2}(\partial_{xq} e_i^p e_i^q - e_i^p \partial_{xq} e_i^q) + \frac{1}{2} \partial_{xq} (e_i^p e_i^q) = e_i^q \partial_{xq} e_i^p. \quad (3.22) $$

A remapping is thus obtained to write

$$ T_t(x) = x + e_i^q \partial_{xq} e_i \Delta t + e_i \Delta \eta_i = x - \mathbf{w}_S^c \Delta t + (-\mathbf{w}_S^c \Delta t - \boldsymbol{\sigma} \Delta \mathbf{B}), \quad (3.23) $$

already derived for differential 0-form in LU framework (Eq. (3.9)). Therefore, the proposed perturbation mapping can also encompass the LU framework for n-forms, and its capacity—given by the Reynolds transport theorem—to deal with extensive properties.

Moreover, for incompressible flows, LU equation further imposes that

$$ \left\{ \begin{array}{l} \nabla \cdot \boldsymbol{\sigma} = 0, \\ \nabla \cdot \nabla \cdot \mathbf{a} = 0. \end{array} \right. \quad (3.24) $$

Translating it into our present notation, it reads as

$$ \left\{ \begin{array}{l} \partial_{x_p} e_i^p = 0, \text{ for each } i \\ \partial_{x_p} \partial_{x_q} (e_i^p e_i^q) = 0. \end{array} \right. \quad (3.25) $$

Following straightforward calculation, Eq. (3.24) is found equivalent to that $T_t^* \theta = \theta$ for $\theta = dx^1 \wedge \cdots \wedge dx^n$. Such a result is expected since constraints Eq. (3.24) are obtained from the LU density conservation.
## 3.2 The SALT Perturbation Scheme

Holm (2015) derived the original SALT equation following a stochastically constrained variational principle $\delta S = 0$, for which

$$ \left\{ \begin{array}{l} S(u, q) = \int \ell(u, q) \partial t, \\ \partial q + \ell_{\partial x_i} q = 0, \end{array} \right. \qquad (3.26) $$

where $\ell(u, q)$ is the Lagrangian of the system, $\ell$ is the Lie derivative, and $x_t(x)$ is defined by (using our notation)

$$ x_t(x) = x_0(x) + \int_0^t u(x, s) \partial s - \int_0^t e_i(x) \circ \partial \eta_i(s), \qquad (3.27) $$

in which $u$ is the velocity vector field. The $\circ$ means that the integral is defined in the Stratonovich sense, instead of in the Ito sense. Hence, $\partial x_t = u(x, t) \partial t - e_i \circ \partial \eta_i$ refers to an infinitesimal stochastic tangent field on the domain. We can express $\partial x_t = T_t(x) - x + u \partial t$. Note the difference between Ito's notation and Stratonovich's notation, i.e. $e_i \circ \partial \eta_i \neq e_i \partial \eta_i$. The initial expression of $T_t$ essentially follows Ito's notation. In this subsection, it comes that $T_t(x) \neq x - e_i \Delta \eta_i$. Instead, it becomes $T_t(x) = x + \frac{1}{2} e_i^p \partial_{x_p} e_i \Delta t - e_i \Delta \eta_i$.

In the second equation of Eq. (3.26), $q$ is assumed to be a quantity advected by the flow. $q$ can correspond to any differential form that is not uniquely determined by the velocity (since the SALT equation for the velocity is usually determined by the first equation of Eq. (3.26)). Holm (2015) evaluates the Lie derivative $\ell_{\partial x_i} q$ using Cartan's formula:

$$ \ell_{\partial x_i} q = d(i_{\partial x_i} q) + i_{\partial x_i} dq. \qquad (3.28) $$

This Lie derivative $\ell_{\partial x_i} q$ corresponds to $T_t^* q - q + f^q(S) \partial r$, if we assume that the deterministic forecast of $q$ is simply the advection of $q$ by $u$. More generally, $\ell_{\partial x_t - u \partial r} q = T_t^* q - q$. Therefore, the SALT equation for $q$ is the same as our perturbation for $q$. Note, the Cartan's formula can not be directly applied to calculate the Lie derivative if the expression of $\partial x_t$ is in Ito's notation.

Within the SALT setting, the velocity $u$ comes from the first equation of Eq. (3.26). For most cases, the velocity $u$ is associated with the momentum, a differential 1-form $m = u^j dx^j = u^1 dx^1 + \dots + u^n dx^n$. When the Lagrangian includes the kinetic energy, Holm (2015) observed that the stochastic noises contribute a term $\ell_{\partial x_i} \theta$, where $\theta$ is a differential 1-form related to the momentum 1-form. In particular, $\theta = m$ for the “Stratonovich stochastic Euler-Poincaré flow” example, and $\theta = m + R^j dx^j$ for the “Stochastic Euler-Boussinesq equations of a rotating stratified incompressible fluid”.
Already pointed out, the operator $\mathcal{L}_{\partial x_i}$ is closely related to $T_i^*$, and the SALT momentum equation can thus also be derived using our proposed perturbation scheme by properly choosing $\theta$, without relying on Lagrangian mechanics.

Another way to appreciate the correspondence to SALT is by looking at the final SPDE. If we choose $\theta$ to be a differential 1-form to represent the momentum $f$, i.e. $\theta = f^j dx^j$ we obtain from Eq. (1.7) (more details in Zhen et al 2022):

$$ \begin{aligned} \mathcal{L}_s f^j = & (a^p \partial_{x^p} f^j + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f^j + \partial_{x^j} a^p f^p + \partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p) \mathcal{L}t \\ & + (e_i^p \partial_{x^p} f^j + \partial_{x^j} e_i^p f^p) \mathcal{L}\eta_i. \end{aligned} \quad (3.29) $$

Regrouping the terms for physical interpretation, it writes:

$$ \begin{aligned} \frac{\partial_s f^j}{\partial t} = & -V^p \partial_{x^p} f^j + \partial_{x^p} \left( \left( \frac{1}{2} e_i^p e_i^q \right) \partial_{x^q} f^j \right) \\ & + \partial_{x^j} \left( a^p + e_i^p \frac{\partial \eta_i}{\partial t} \right) f^p + \partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p. \end{aligned} \quad (3.30) $$

Two last terms of the right-hand side complete the advection-diffusion terms, already appearing in (3.7). The first one, $\partial_{x^j} \left( -a^p - e_i^p \frac{\partial \eta_i}{\partial t} \right) f^p$, is reminiscent to the additional terms appearing in SALT momentum equations (Holm 2015, Resseguier et al 2020b). The second term, $-\partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p$, comes from cross-correlation in Itô notation.

# 4 Conclusion

As demonstrated, both SALT and LU equations can be recovered using a prescribed definition of a random diffeomorphism $T_t$ used to perturb the physical space. However, compared with SALT and LU settings, the proposed perturbation scheme does not directly rely on a particular physics. Hence, the random mapping is more flexible and can be applied to any PDE. Interestingly, similarities and differences can then be identified and studied between the proposed use of the random diffeomorphism and the existing stochastic physical SALT and LU settings. For instance, the proposed derivation provides an interesting interpretation the operator $\mathcal{L}_{\partial x_t - u \partial t}$, appearing in the SALT equation. This term can indeed represent an infinitesimal forecast error at every forecast time step.

To apply the proposed perturbation scheme to any specific model, the diffeomorphism parameters $a$ and $e_i$ must be determined specifically. Hence it is necessary to learn these parameters from existing data, experimental runs, or additional physical considerations. This framework naturally provides new perspectives to generate ensembles through constrained stochastic mappings applied in the physical space.
**Acknowledgments** This work is supported by the ERC project 856408-STUOD and SCALIAN DS.

# Appendix: Expression of $T_t^*\theta$

Given coordinates $(x^1, \dots, x^n)$ and a differential $k$-form $\theta$, Zhen et al (2022) (Appendix B) proof that:

$$ 
\begin{aligned}
T_t^*\theta &= f(T_t(x))T_t^*(dx^{i_1} \wedge \dots \wedge dx^{i_k}) \\
&= \theta + \left\{ \left(\langle \nabla f, a \rangle + \frac{1}{2} e_i^\top H_f e_i + I W\right) dx^{i_1} \wedge \dots \wedge dx^{i_n} \right. \\
&\qquad + \sum_{s=1}^k f \partial_{x^j} a^{i_s} dx^{i_1} \wedge \dots dx^j \wedge \dots \wedge dx^{i_k} \\
&\qquad + \left( \sum_{s<r} f \partial_{x^j} e_i^{i_s} \partial_{x^l} e_i^{i_r} dx^{i_1} \wedge \dots \wedge dx^j \wedge \dots \wedge dx^l \wedge \dots \wedge dx^{i_k} \right) \\
&\qquad + \left( \sum_{s=1}^k \langle \nabla f, e_i \rangle \partial_{x^j} e_i^{i_s} dx^{i_1} \wedge \dots \wedge dx^j \wedge \dots \wedge dx^{i_k} \right) \right\} \Delta t \\
&\qquad + \left\{ \langle \nabla f, e_i \rangle dx^{i_1} \wedge \dots \wedge dx^{i_k} \right. \\
&\qquad \left. + \sum_{s=1}^k f \partial_{x^j} e_i^{i_s} dx^{i_1} \wedge \dots \wedge dx^j \wedge \dots \wedge dx^{i_k} \right\} \Delta \eta_i \\
&\qquad + o(\Delta t),
\end{aligned} 
\quad (A.1) 
$$

where $IW$ is the additional term appearing in the Itô-Wentzell formula (Kunita 1997). Here, there is no noise in the original dynamics (2.3) (the first step (2.4) of the randomized dynamics) which could be correlated with the noise of the resulting stochastic scheme (2.5). That is why $IW = 0$ in the above Taylor development of $f$. Indeed, there is no additional cross-correlation term between $T_t^*$ and $\tilde{\theta}(t + \Delta t) = \theta(t) + g^0(S(t))\Delta t$. The final SPDE (2.8) makes clear the link between the solution $\theta$ and the Brownian motions $\eta_i$. But, at a given time step $t$, since (2.2) has no noise term, $\tilde{\theta}(t + \Delta t)$ is correlated with the $t' \mapsto \eta_i(t')$ for $t' < t$ only, and is independent of the new Brownian increment $\Delta \eta_i(t)$ generating $T_t$. Therefore, there is no cross-correlation term between $T_t^*$ and $\tilde{\theta}(t + \Delta t)$.

To simplify Eq. (A.1), wedge algebra is applied and the high-order infinitesimal $o(\Delta t)$ is ignored. Accordingly, $T_t^*\theta$ is more compactly written as

$$ T_t^*\theta = \theta + \mathcal{M}(\theta)\Delta t + \mathcal{N}_i(\theta)\Delta \eta_i, \quad (A.2) $$

for some differential $k$-forms $\mathcal{M}(\theta)$ and $\mathcal{N}_i(\theta)$.
## References

1. Bethencourt De Leon A (2021) On the effect of stochastic lie transport noise on fluid dynamic equations. PhD thesis, Imperial College London

2. Dufée B, Mémin E, Crisan D (2022) Stochastic parametrization: an alternative to inflation in ensemble kalman filters. Quarterly Journal of the Royal Meteorological Society 148(744):1075–1091

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

4. Kunita H (1997) Stochastic flows and stochastic differential equations, vol 24. Cambridge university press

5. Mémin E (2014) Fluid flow dynamics under location uncertainty. Geophysical & Astrophysical Fluid Dynamics 108(2):119–146. https://doi.org/10.1080/03091929.2013.836190

6. Mikulevicius R, Rozovskii B (2004) Stochastic Navier–Stokes equations for turbulent flows. SIAM Journal on Mathematical Analysis 35(5):1250–1310

7. Resseguier V (2017) Mixing and fluid dynamics under location uncertainty. PhD thesis, Rennes 1 university

8. Resseguier V, Mémin E, Chapron B (2017) Geophysical flows under location uncertainty, part I random transport and general models. Geophysical & Astrophysical Fluid Dynamics 111(3):149–176

9. Resseguier V, Li L, Jouan G, Dérian P, Mémin E, Bertrand C (2020a) New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics. Archives of Computational Methods in Engineering pp 1–82

10. Resseguier V, Pan W, Fox-Kemper B (2020b) Data-driven versus self-similar parameterizations for stochastic advection by lie transport and location uncertainty. Nonlinear Processes in Geophysics 27(2):209–234

11. Tandeo P, Ailliot P, Bocquet M, Carrassi A, Carrassi A, Miyoshi T, Pulido M, Pulido M, Zhen Y (2020) A review of innovation-based methods to jointly estimate model and observation error covariance matrices in ensemble data assimilation. Monthly Weather Review

12. Zhen Y, Resseguier V, Chapron B (2022) Physically constrained covariance inflation from location uncertainty and optimal transportation. arXiv preprint arXiv:221104207

**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_116.png" /></figure>
<figure><img src="image_117.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 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( 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_σ \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_σ \cdot dB_t) \, dx. \quad (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_σ \cdot dB_t. \quad (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_σ) dt - \frac{1}{2} \nabla \cdot (a \nabla q) dt = Q_t dt + Q_σ \cdot dB_t. \quad (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( 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 (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 (\mathbf{a} \nabla (\rho u_i)) 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 \\ & \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
$$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}} \\ & + \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} \quad (23) $$
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}
& \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} = \mathcal{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
$$ \rho = \rho_0 \left(1 - \beta_T (T - T_0) + \beta_p p\right), \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
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) \mathrm{d}z \tag{47} \\ 
d p_t^\sigma &= -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \mathrm{d}z. \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} \mathrm{d}t - \frac{1}{\rho_0} \frac{\partial \mathrm{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) \mathrm{d}z \\
\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) \mathrm{d}z \\
d p_t^\sigma &= -\rho_0 \int_z^\eta (\sigma_t d\mathbf{B}_t \cdot \nabla) w \mathrm{d}z \\
\mathbb{D}_t b + (w^* \mathrm{d}t + (\sigma_t d\mathbf{B}_t)_z) N^2 - \frac{1}{2} \nabla \cdot (a_\bullet z N^2) \mathrm{d}t &= \frac{g\beta_T}{\gamma\rho_0} (dQ + dW_d) \\
dQ &= dQ_{\text{rad}} + \frac{1}{RePr} \nabla \cdot (\nabla T) \mathrm{d}t \\
dW_d &= \left( \frac{1}{2} \nabla \cdot \mathbf{a} \mathrm{d}t - \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}{\partial x_i} p_t^\sigma && \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 \\
dp_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 (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.
# 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) $$
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
$$ \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. $$