thinkings on the stochastic noise. The left-hand side terms corresponds to the
vertical acceleration with a scaling such that the noise vertical acceleration is
strong compared to the large-scale vertical acceleration terms. Note that the two
underlined terms in (7) are Brownian terms emerging from the stochastic pressure
formulation (6) on the right-hand side and from the vertical velocity transport by
the noise on the left-hand side. The two other blue terms on the LHS are associated
with modified advection and dissipation (projected on the vertical velocity $w$) that
emerged through the three dimensional generalization of (4):

$$
\mathbb{D}_t \mathbf{u} = d_t \mathbf{u}_h + \mathbf{u} \cdot \nabla \mathbf{u}_h + \left( -\frac{1}{2} \nabla \cdot \mathbf{a} dt + \sigma d\mathbf{B}_t \right) \cdot \nabla \mathbf{u} - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \mathbf{u}) dt, \quad (8)
$$

where the material derivative of vertical velocities associated with the resolved flow
(i.e. $d_t w + \mathbf{u} \cdot \nabla w$) has been neglected.

The noise being given (and calibrated from data or a known relation), from (7), it
is thus possible to compute the various Brownian terms on the LHS, then to integrate
vertically the results to obtain a 3D map of the modified pressure field as a result of
the noise transport. Separating safely the martingale part (Brownian terms) from the
smooth-in-time components (“dt” terms), we have

$$
d p_t^\sigma(z) = d p_t^\sigma|_{z=\eta} + \rho_0 \int_z^\eta (\sigma d\mathbf{B}_t \cdot \nabla w) dz', \quad (9)
$$

for the martingale component, and

$$
p(z)dt = p|_{z=\eta}dt + \rho_0 g(\eta - z)dt \\
+ \rho_0 \int_z^\eta \left( bdt - \underbrace{\left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \right) \cdot \nabla w dt + \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla w) dt \right)}_{b_{NH}} \right) dz', \quad (10)
$$

for the smooth-in-time component. The three last terms on the RHS of (10) can be
compared to the deterministic non-hydrostatic pressure correction of Klingbeil and
Burchard (2013), although the material derivative of $w$ associated with resolved flow
is not included in our stochastic formulation. It can be noted that our formulation
involves a 3D diffusion of the vertical velocity ensuing from the noise action
as well as the contribution of the modified Ito-Stokes term arising from the
spatial inhomogeneity of the noise. Both (9) and (10) should be integrated with
appropriate boundary conditions at $\eta$ to incorporate fast and smooth-in-time surface
pressure contributions (such as surface waves or atmospheric pressure loading,
respectively), but such contributions can be neglected at first approximation without
loss of generality. Results of (9) and (10) can then be used to feedback onto the
horizontal momentum equation (3) solved by an hydrostatic model. Assuming a
strict separation of the martingale part and the smooth-in-time component, only (10)
is assumed to feedback onto the resolved flow. The martingale components are assumed to balance each other, thus not affecting the resolved flow. This assumption can be interpreted as a Large Eddy Simulation (LES) -like approach, as discussed by Bauer et al. (2020).

As a preliminary step, we will further simplify the structure of the variance tensor $a$ in order to reduce the second and third terms on the RHS of (10) to a simple Laplacian viscosity—induced here by the noise contribution. This simplification is motivated in the following. In the LU framework, the strength of the noise is measured by its (one-point co-) variance, such that

$$a(x, t) \stackrel{\Delta}{=} \check{q}(x, x, t), \qquad (11)$$

with $\check{q}(x, x, t)$ a matrix kernel defined as

$$\check{q}(x, x, t) \stackrel{\Delta}{=} \int_{\Omega} \check{\sigma}(x, x', t) \check{\sigma}(x, x', t)^{T} dx', \qquad (12)$$

with $\check{\sigma}(\cdot, \cdot, t)$ a bounded matrix kernel defining the correlation deterministic integral operator $\sigma_t : L^2(\Omega) \to L^2(\Omega)$

$$\sigma_t f(x) \stackrel{\Delta}{=} \int_{\Omega} \check{\sigma}(x, y, t) f(y) dy, \quad \forall f \in (L^2(\Omega)). \qquad (13)$$

(See, e.g. Bauer et al. 2020, Mémin 2014, Resseguier et al. 2017, for further details). Although the previous definition of the noise is general, it is possible, through the Mercer's theorem, to express the noise variance as a spectral decomposition of the form:

$$a(x, t) = \sum_{n \in \mathbb{N}} \lambda_n(t) \phi_n(x, t) \phi_n^T(x, t), \qquad (14)$$

where $\phi_n(x, t)$ define an orthonormal eigenfunction basis of the correlation operator, $\sigma_t$, with $\lambda_n(t)$ their corresponding eigenvalues. For a stationary noise, this reduces to a classical POD (or EOF) decomposition, in which the eigenfunctions are the solution of the eigenvalue problem

$$\int_{\Omega} K(x, x') \phi_n(x') dx' = \lambda_n \phi_n(x) \qquad (15)$$

with $K$ the two-point correlation tensor.

The next step is to assume isotropy and homogeneity of the noise structure, in which case the Fourier modes $\phi_n = e^{2\pi i k \cdot x}$ are a natural choice to satisfy (15), which implies (Berkooz et al. 1993):

$$K = \sum_{n} \lambda_{n} e^{2\pi i k \cdot x} e^{-2\pi i k \cdot x'}. \qquad (16)$$
Under isotropic condition, the variance of the divergence-free noise is constant and diagonal, such that the first term associated with $b_{NH}$ in (10) is identically zero, and the noise induced dissipation reduces to:

$$ \frac{1}{2} \nabla \cdot (a \nabla w) = \nu \Delta w, \qquad (17) $$

with $\nu$ the isotropic, homogeneous noise induced momentum dissipation. Through vertical integration of (17), we recover part of the non-hydrostatic pressure correction proposed by Klingbeil and Burchard (2013), which in the present case mimic entrainment/detrainment of convective plumes leading to changes in their spatial organization. This modified HPE will be termed Quasi-Nonhydrostatic (Q-NH), by analogy with the Quasi-Hydrostatic (QH) of Marshall et al. (1997) where non-traditional Coriolis terms are added into the HPE.

# 3 Numerical Implementation and Simulations

Our objective is to implement this stochastic, non-hydrostatic pressure correction in the hydrostatic kernel of the Coastal and Regional Ocean Community model (CROCO; http://www.croco-ocean.org). CROCO is a new ocean model that builds upon the structure of the ROMS-AGRIF primitive equation solver (Shchepetkin and McWilliams 2005, Debreu et al. 2012). The non-hydrostatic, non-Boussinesq (NQB; Auclair et al. 2018) capabilities of CROCO will also be used to construct a reference simulation for validation (see Table 1). We review in the following some important steps for the implementation of the stochastic pressure correction within CROCO, discuss their implications and how we treat the pressure correction within the hydrostatic CROCO kernel.

## 3.1 Stochastic, Non-hydrostatic Pressure Correction

In its hydrostatic mode, CROCO computes $\frac{p}{\rho_0}$, from which horizontal gradients directly feed the baroclinic horizontal momentum equation (i.e. (2a) and (2b)). Our strategy is to include the NH pressure correction via a modified density/buoyancy field, such that the pressure field becomes:

$$ \frac{p(z)}{\rho_0} = \frac{p|_{z=\eta}}{\rho_0} + \int_{z}^{\eta} \frac{(\rho - \rho_0) - b_{NH}}{\rho_0} g \, dz', \qquad (18) $$

with $b_{NH}$ collecting the different contributions of the vertical momentum equation ($w_trends$) contributing in the pressure correction, normalized by gravity:
$$b_{NH} = \frac{1}{g} \sum w_{\text{trends}}, \qquad (19)$$

where the horizontal dissipation of vertical velocity (e.g. Eq. (17)) is computed along sigma coordinates. Our strategy is similar to Delorme et al. (2021), who derived a Quasi-Hydrostatic version of CROCO by including the non-traditional Coriolis effects through buoyancy correction.

Let us note that $b_{NH}$ is abusively denoted through a buoyancy variable, however it corresponds to corrections brought by the noise to the usual hydrostatic pressure. Such a correction should not be interpreted as an actual modification of the density/buoyancy field of the stratified ocean. Thus, the stochastic contribution is not included in the specific treatment of baroclinic-barotropic mode coupling of CROCO, which aims at accounting for the non-uniform density field for the propagation of gravity waves (Gill 1982), ultimately reducing the usual mode-coupling error associated with mode-splitting schemes (Shchepetkin and McWilliams 2005). In other words, we do not expect this ‘non-hydrostatic buoyancy’ to affect gravity wave’s propagation.

Finally, CROCO uses a third-order predictor-corrector (LF-AM3) time-stepping scheme for tracers and baroclinic momentum. This scheme consists of a Leapfrog (LF) predictor with 3rd-order Adams-Moulton (AM) interpolation. It also uses split-explicit techniques to robustly couple the slow, baroclinic and the fast, barotropic modes associated with the time evolving non-linear free surface. A complete description of the several stages of CROCO time-stepping can be found in Section 5 of Shchepetkin and McWilliams (2005). This predictor-corrector, split-explicit scheme implies some tendency terms of the baroclinic mode are computed twice to step forward baroclinic momentum and tracer equations from time step $t$ to time step $t + \Delta t$. The first computation is performed at the prediction stage, and the second computation is performed at the correction stage. These tendencies include pressure gradients. To avoid double counting the stochastic pressure correction and for stability reasons, the modified non-hydrostatic buoyancy is computed only at the correction stage.

## 3.2 Numerical Experiments

The numerical experiments we used to test our stochastic non-hydrostatic pressure correction are oceanic deep convection events. The configuration is inspired by free convection studies (e.g. Souza et al. 2020) where an horizontally uniform surface cooling is applied to a constantly stratified, horizontally uniform ocean. Although our interest is on deep ocean convection rather than mixed layer free convection, we have adopted an horizontally uniform setting (as usually done in free convection) instead of a horizontally structured system as proposed earlier by, e.g. Marshall and Schott (1999), where surface cooling is confined within a specified region (usually a disc). This setting allows the analysis of interacting
<figure><img src="image_41.png" /><figcaption>**Fig. 1** Illustration of the 3D structure of the simulation after three days of simulation in a non-hydrostatic setting. The color shading is temperature</figcaption></figure>

convective plumes with non-convective environment. However, such configurations are usually conducted at coarser resolution and oriented toward process studies of the geostrophic organization of convective plumes. Here, our focus is on parameterization and we wish to start with simplified settings in order to capture the essence of deep convection dynamics; interactions with a prescribed background, non-convective environment is left for further work.

From this horizontally uniform and vertically constant stratified initial condition, the model is stepped forward in time on a 100×100×100 discretized grid points with isotropic resolution of 10 meters, and exposed to a constant (in time and space) cooling rate of $Q_T = -500 \text{ W m}^{-2}$ heat flux (Fig. 1). This leads to a cooling of upper ocean layers, which ultimately become unstable through static instabilities as a result of a negative buoyancy frequency ($N^2 = -\frac{g}{\rho}\partial_z b) < 0$, thus undergoing convection. The current settings are run with no Coriolis forcing, i.e. $f = 0 \text{ s}^{-1}$; inclusion of Coriolis effects will be the subject of further work. The model is initialized with stochastic perturbations on ocean upper layers temperature decaying with depth (following Souza et al. 2020) to trigger the formation of convective plumes

$$T(x, y, z)|_{t=0} = T(z) + \sum_{(m,n)=0}^{10} \left( e^{2\pi(\mathbf{k}\cdot\mathbf{x}+\phi_{n,m})} \right) \mathcal{N}(0, 1) * \sqrt{\sigma^2} * e^{40z/Nz} \quad (20)$$

with $T(z) = T|_{z=0} - \alpha z$, $T|_{z=0} = 3 \text{ K}$ and $\alpha$ is a constant defined as $\alpha = \frac{1.9e^{-6}}{g*(\alpha_T/\rho_0)}$ ($\alpha_T = 0.2048 \text{ K}^{-1}$ is the thermal coefficient, $g = 9.81 \text{ m s}^{-2}$ is gravity and $\rho_0 = 1024 \text{ kg m}^{-3}$ is reference density). The second term on the RHS of (20) is the stochastic perturbation defined as the sum of plane waves with random phase $\phi_{m,n}$
<table><caption>**Table 1** Summary of the experiments and their numerical details. NBQ stands for the non-hydrostatic, non-Boussinesq CROCO kernel of Auclair et al. (2018); NH, Hydro and Q-NH stand for Non-Hydrostatic, Hydrostatic and Quasi-Nonhydrostatic; WENO5 and C4 for the 5-th order and the 4th-order centred advection schemes; KPP for the K-Profil Parameterization of Large et al. (1994)</caption><thead><tr><th>Name</th><th>Kernel</th><th>Hz adv</th><th>Vert. adv</th><th>Closure</th><th>(Δx, Δz)</th><th>Δt</th><th>ν<sub>w</sub></th></tr></thead><tbody><tr><td>NH</td><td>NBQ</td><td>WENO5</td><td>WENO5</td><td>–</td><td>(10, 10 m)</td><td>2.5 s</td><td>(implicit)</td></tr><tr><td>Hydro</td><td>Hydro</td><td>WENO5</td><td>C4</td><td>KPP</td><td>(10, 10 m)</td><td>2.5 s</td><td>–</td></tr><tr><td>Q-NH</td><td>Hydro</td><td>WENO5</td><td>C4</td><td>KPP</td><td>(10, 10 m)</td><td>2.5 s</td><td>1 m<sup>2</sup>s<sup>-1</sup></td></tr></tbody></table>

and of amplitude $N * \sqrt{\sigma^2}$, where $N$ is a Gaussian white noise distribution and $\sigma^2 = 10^{-8} K^2$ represents the variance of the stochastic perturbations. The random phases are drawn from an uniform distribution over the range [0, 1].

This configuration has been integrated forward in time to produce several numerical experiments in order to assess the performance of our implementation. It includes a pure Non-Hydrostatic (NH), which make use of the non-hydrostatic non-Boussinesq capabilities of CROCO (NBQ, Auclair et al. 2018), and a pure hydrostatic (Hydro) reference experiments. We then compare the solution produced by our Quasi-Nonhydrostatic (Q-NH) experiment, which includes the stochastic pressure correction, with the solutions produced by Hydro and NH.

For both Q-NH and Hydro, the KPP (Large et al. 1994) closure scheme is used to represent vertical sub-grid scale mixing. As stated in introduction, this scheme mimics vertical fluxes through a dissipative down-gradient operator. In the case of static instability associated with convective events, the dissipation coefficient is set to $K_\theta = 0.1 \text{ m}^2\text{s}^{-1}$ in CROCO. Sensitivity tests (not shown) using TKE (Gaspar et al. 1990) closure scheme instead revealed that the choice of the closure scheme has little effect on the solution produced by our Hydro experiment.

Horizontal dissipation of vertical velocity in Q-NH is set to $v_w = 1 \text{ m}^2\text{s}^{-1}$, which corresponds to high values of dissipation estimated through a Smagorinsky-like approach $v_{Smago} = \frac{\alpha}{2} \Delta_{xy}^2 ((\partial_x w)^2 + (\partial_y w)^2)$, with $\Delta_{xy} = 10 \text{ m}$ the horizontal resolution of our configuration, and $\alpha = 0.2$. Table 1 summarizes the different experiments, along with some numerical details.

Finally, note that all three experiments are conducted at the same isotropic resolution of 10m. Evaluating the performance of our Q-NH model for climate scale regimes (i.e. with horizontal resolution much coarser than vertical resolution) will be the subject of further work.

# 4 Results

We show on Fig. 2 snapshots of the vertical velocities as simulated by NH, Hydro and Q-NH after 3 days of simulation. Obviously, the NH experiment produces weaker and larger scale structures as compared to the two other experiments. It
<figure><img src="image_42.png" /><figcaption>**Fig. 2** Horizontal (top) and vertical (bottom) sections of vertical velocities for the Non-Hydrostatic (NH, left), the Hydrostatic (Hydro, center) and the Quasi-Nonhydrostatic (Q-NH, right) run. Snapshots are shown after 3 days of simulation with all other components (forcing, dissipation, stratification, resolution) held constant</figcaption></figure>

is noticeable, however, that the amplitude of *w* in Q-NH is reduced as compared to Hydro, with larger scale structures of the convective plumes. In the Hydro experiment, plumes are localized near the grid scale and exhibit almost an order of magnitude larger vertical velocities as compared to the NH reference. The reduced vertical velocities in Q-NH and broadening of the associate spatial scale of the plumes can be interpreted as a result of the entrainment/detrainment mechanism, which is here simply represented as a purely horizontal viscous stress on vertical velocities.

Aside from vertical velocities, it is also instructive to analyse the consequence for the temperature profile, an indication of the capability of convection in producing deep water masses. Figure 3 shows the horizontally averaged temperature vertical profile, along with the rate of the Mixed Layer Depth (MLD) deepening, for the different experiments. A first remarkable result is the similarity between the temperature profiles produced by all the experiments within the MLD (i.e. $z < -300$ m). Additional tests (not shown) indeed reveal the very weak sensitivity to numerical implementation (i.e. non-hydrostatic, hydrostatic with different vertical mixing schemes), as well as the level of dissipation. The most noticeable differences between the experiments appear at the base of the mixed layer. In particular, the consequences of too strong vertical velocities in Hydro is to produce too deep water masses with a too strong penetrative convection (i.e. the envelop at the base of the mixed layer where water masses are warmer than their initial state). Although vertical velocities in Q-NH remain significantly larger than those produced by NH, the effects of horizontal viscous forces on the vertical velocities is to significantly damp the penetration of convective plumes below the mixed layer, inducing a strong reduction of water masses formation. The reduction of penetrative convection induces significant biases in the deepening rate of the MLD (Fig. 3, right panel).
<figure><img src="image_43.png" /><figcaption>**Fig. 3** Horizontally averaged temperature profiles (left) and deepening rate of the mixed layer depth (MLD, right) for the different runs at the end of the 3-day long simulations. Deepening MLD rate are compared to the analytical estimates of Marshall and Schott (1999) and Souza et al. (2020)</figcaption></figure>

In Q-NH, MLD seats between the theoretical predictions of Marshall and Schott (1999) and that of Souza et al. (2020), where the former do not consider penetrative convection in their scaling while the latter do. That MLD is too shallow in Q-NH, as compared to NH, is likely a consequence of too strong dissipation imposed to the system. We note, however, that we have only considered the dissipative, rectification contribution (i.e. smooth-in-time) of the stochastic transport as a result of the strict separation assumption between martingale and smooth-in-time components. Further work are required to evaluate how the Brownian part of the stochastic transport impact the deepening of the MLD.

# 5 Conclusion and Perspectives

In this study, we detailed the first steps toward a full stochastic parameterization of deep ocean convection along with their implementation in the general circulation model CROCO. Our preliminary results, which consist of an approximation of the horizontal noise structure as homogeneous and isotropic, led us to recover part of the derivation provided by Klingbeil and Burchard (2013) in a deterministic case. Our results are encouraging, and we are now in a position of extending the current analysis to a fully consistent stochastic framework.

The first step in this direction will be to implement the stochastic pressure noise contribution which comes in pair with the idealized Laplacian horizontal viscosity action on vertical velocities in the context of an hydrostatic simulation of deep ocean convection (i.e. Eq. (9)). Following previous work of Pierre Dérian and
Etienne Mémin ('Hyper-viscosity' noise for transport under location uncertainty), we will consider a simplified expression of the stochastic transport of vertical velocity to construct our stochastic pressure noise. This approach is meant to obtain the stochastic transport associated with Laplacian or hyper-viscosity dissipation as usually implemented in OGCM. With these considerations, it is possible to express the stochastic transport of (9) as:

$$ \sigma d\mathbf{B}_t \cdot \nabla w = \sum_k \gamma_k \lambda_k e_k(\mathbf{x}) \quad (21) $$

with $e_k$ a basis, defined here as Daubechies wavelets, $\gamma_k$ denotes independent normally distributed variables, and $\lambda_k$ are the wavelet coefficients defined as:

$$ \lambda_k = \left\langle \sqrt{2\epsilon dt} | v^{1/2} \nabla w |; e_k \right\rangle_{L^2}, \quad (22) $$

with $\epsilon$ a scaling factor controlling the ratio of variance created by the noise to energy dissipation. Accounting for the stochastic pressure would also require considerations for the Brownian components of the stochastic transport in the horizontal momentum advection. These steps are part of further works to achieve a full, consistent implementation of LU transport in CROCO.

**Acknowledgments** This work is supported by the ERC project 856408-STUOD.

## References

* Auclair F, Bordois L, Dossmann Y, Duhaut T, Paci A, Ulses C, Nguyen C (2018) A non-hydrostatic non-Boussinesq algorithm for free-surface ocean modelling. Ocean Model 132:12–29

* Bauer W, Chandramouli P, Chapron B, Li L, Mémin E (2020) Deciphering the role of small-scale inhomogeneity on geophysical flow structuration: a stochastic approach. J Phys Oceanogr 50(4):983–1003

* Berkooz G, Holmes P, Lumley JL (1993) The proper orthogonal decomposition in the analysis of turbulent flows. Annual review of fluid mechanics 25(1):539–575

* Debreu L, Marchesiello P, Penven P, Cambon G (2012) Two-way nesting in split-explicit ocean models: Algorithms, implementation and validation. Ocean Model 49:1–21

* Delorme BL, Thomas LN, Marchesiello P, Gula J, Roullet G, Molemaker MJ (2021) Enhanced abyssal mixing in the equatorial pacific associated with non-traditional effects. J Phys Oceanogr 51(6):1895–1914

* Gaspar P, Grégoris Y, Lefevre JM (1990) A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station papa and long-term upper ocean study site. Journal of Geophysical Research: Oceans 95(C9):16,179–16,193

* Gill AE (1982) Atmosphere-ocean dynamics, vol 30. Academic press

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

* Heuzé C (2017) North atlantic deep water formation and amoc in cmip5 models. Ocean Sci 13(4):609–622

* Heuzé C (2021) Antarctic bottom water and north atlantic deep water in cmip6 models. Ocean Sci 17(1):59–90
Hourdin F, Musat I, Bony S, Braconnot P, Codron F, Dufresne JL, Fairhead L, Filiberti MA, Friedlingstein P, Grandpeix JY, et al (2006) The Imdz4 general circulation model: climate performance and sensitivity to parametrized physics with emphasis on tropical convection. Clim Dyn 27:787-813

Klingbeil K, Burchard H (2013) Implementation of a direct nonhydrostatic pressure gradient discretisation into a layered ocean model. Ocean Model 65:64-77

Large WG, McWilliams JC, Doney SC (1994) Oceanic vertical mixing: A review and a model with a nonlocal boundary layer parameterization. Rev Geophys 32(4):363-403

Marshall J, Schott F (1999) Open-ocean convection: Observations, theory, and models. Reviews of Geophysics 37(1):1-64

Marshall J, Hill C, Perelman L, Adcroft A (1997) Hydrostatic, quasi-hydrostatic, and nonhydrostatic ocean modeling. Journal of Geophysical Research: Oceans 102(C3):5733-5752

Masson-Delmotte V, Zhai P, Pirani A, Connors SL, Péan C, Berger S, Caud N, Chen Y, Goldfarb L, Gomis M, et al (2021) Climate change 2021: the physical science basis. Contribution of working group I to the sixth assessment report of the intergovernmental panel on climate change 2

Mémin E (2014) Fluid flow dynamics under location uncertainty. Geophys Astrophys Fluid Dynamics 108(2):119-146

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

Shchepetkin AF, McWilliams JC (2005) The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean modelling 9(4):347-404

Souza AN, Wagner G, Ramadhan A, Allen B, Churavy V, Schloss J, Campin J, Hill C, Edelman A, Marshall J, et al (2020) Uncertainty Quantification of Ocean Parameterizations: Application to the K-Profile-Parameterization for Penetrative Convection. Journal of Advances in Modeling Earth Systems 12(12):e2020MS002,108

Suselj K, Kurowski MJ, Teixeira J (2019) A unified eddy-diffusivity/mass-flux approach for modeling atmospheric convection. J Atmos Sci 76(8):2505-2537

Tucciarone FL, Mémin E, Li L (2023) Primitive equations under location uncertainty: Analytical description and model development. In: Stochastic Transport in Upper Ocean Dynamics Annual Workshop, Springer, pp 287-300

Umlauf L, Burchard H (2003) A generic length-scale equation for geophysical turbulence models. J Mar Res 61(2):235-265

Zhang R, Sutton R, Danabasoglu G, Kwon YO, Marsh R, Yeager SG, Amrhein DE, Little CM (2019) A review of the role of the Atlantic meridional overturning circulation in Atlantic multidecadal variability and associated climate impacts. Rev Geophys 57(2):316-375

**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_44.png" /></figure>
# Comparison of Stochastic Parametrization Schemes Using Data Assimilation on Triad Models

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

Check for
updates

Alexander Lobbe, Dan Crisan, Darryl Holm, Etienne Mémin, Oana Lang, and Bertrand Chapron

**Abstract** In recent years, stochastic parametrizations have been ubiquitous in modelling uncertainty in fluid dynamics models. One source of model uncertainty comes from the coarse graining of the fine-scale data and is in common usage in computational simulations at coarser scales. In this paper, we look at two such stochastic parametrizations: the Stochastic Advection by Lie Transport (SALT) parametrization introduced by Holm (Proc A 471(2176):20140963, 19, 2015) and the Location Uncertainty (LU) parametrization introduced by Mémin (Geophys Astrophys Fluid Dyn 108(2):119–146, 2014). Whilst both parametrizations are available for full-scale models, we study their reduced order versions obtained by projecting them on a complex vector Fourier mode triad of eigenfunctions of the curl. Remarkably, these two parametrizations lead to the same reduced order model, which we term the *helicity-preserving stochastic triad* (HST). This reduced order model is then compared with an alternative model which preserves the energy of the system, and which is termed the *energy preserving stochastic triad* (EST). These low-dimensional models are ideal benchmark models for testing new Data Assimilation algorithms: they are easy to implement, exhibit diverse behaviours depending on the choice of the coefficients and come with natural physical properties such as the conservation of energy and helicity.

**Supplementary Information** The online version contains supplementary material available at
https://github.com/lobbe/stochastic-triads.

A. Lobbe (✉) · D. Crisan · D. Holm · O. Lang
Imperial College London, Mathematics, London, UK
e-mail: alex.lobbe@imperial.ac.uk

E. Mémin
Campus Universitaire de Beaulieu, Inria - Institut National de Recherche en Sciences et Technologies du Numérique, Rennes, France

B. Chapron
Ifremer – Institut Français de Recherche pour l’Exploitation de la Mer, Plouzané, France
# 1 Introduction

The introduction of stochasticity in fluid dynamics has recently been the subject of intense research effort. This approach involves using random processes to model, for example, unresolved scales, or to take into account neglected physical effects. A stochastic formulation for the fluid flow introduces a probabilistic basis for modelling unresolved scales. This is different from the deterministic causal modelling which is difficult to achieve in practice due, for instance, to unknown initial conditions. In addition, some phenomena such as energy backscattering are directly accessible as stochastic processes. Another usage of stochastic modelling is to generate ensembles of realizations of the model. This facilitates the analysis of model uncertainty quantification for different low-resolution computational simulations and their usage to approximate the true state of the fluid, instead of using a single high-resolution numerical simulation.

Some stochastic schemes have been proposed in the literature by considering a variety of ad-hoc perturbations. However, a principled approach is desirable. The formulation of stochastic dynamical systems based on physical principles has recently been proposed in various settings. For a review and classification of approaches to stochastic parameterisation based on physical principles, see [4]. The present work treats two additional new approaches. The first one, called stochastic advection by Lie transport (SALT) relies on the variational principle for fluid dynamics [16]. The second one, called modelling under location uncertainty (LU) is derived from Newton's principle [23]. Both frameworks introduce stochasticity into the Lagrangian specification of the flow field, rather than directly into the Eulerian frame.

In the deterministic case, it is known that three-dimensional fluid flows may trigger a cascade of dynamics across multiple length and time scales. This multiscale behavior poses considerable challenges in the computational simulation using standard Navier-Stokes equations (see e.g. [7], [24], [5]). When modelling turbulence numerically, specialised discretisation methods are needed to decompose the underlying partial differential equations into a very large number of ordinary differential equations. Alternative approaches have been introduced where the Navier-Stokes dynamics in the Fourier space is mimicked using a finite number of variables, say $u^1, u^2, \dots, u^N$. The Fourier space is divided into $N$ shells, and each shell $s_i$ comprises the set of wave vectors $s$ with magnitude $|s| \in (s_02^i, s_02^{i+1})$. Each $u^i$ satisfies an ODE and it represents the magnitude of the velocity field on a length scale of $s_i^{-1}$ ([15], [7]). The quadratic nonlinearity in the Navier-Stokes equations produces triads of interacting vector Fourier modes within each shell. Shell models involving multiple triads have had considerable success in modelling energy and helicity cascades, as well as modelling intermittency in chaotic dynamical systems [6, 10, 9]. Simplified shell models with only a few triads date back to the 1970s and have provided major insight into fluid modal interaction. Even the dynamical system representation of Euler's fluid equations on a single triad has been quite insightful, see e.g. [28, 29].
More recently, the problem of correctly parameterizing effects of small-scale physical processes together with the need for probabilistic ensemble forecasting and uncertainty quantification has led to modern stochastic approaches in the study of turbulence using reduced order shell models. In this work we will explore reduced order models for SALT and LU models obtained by projecting onto helical basis functions [6, 28, 29]. These helical basis functions, defined as eigenfunctions of the curl operator, enable one to construct reduced order stochastic models of fluid flow with a simplified nonlinear interaction. As we will see, under projection onto the basis of helical triad modes, both LU and SALT result in the same reduced order model and this projected model conserves helicity, but it does not conserve energy. Because of this coincidence in projecting the SALT and LU models onto the helical basis, a second reduced order scheme with a strong energy conservation property inspired by [17] and known as the energy preserving stochastic triad (EST) model will be proposed for comparison.

While the EST model is not of transport type, it will provide comparison between two different classes of stochastic dynamical systems. The two classes treated here are (1) the helicity preserving stochastic triad (HST) (comprising both LU and SALT on the helical basis) and (2) the energy preserving stochastic triad (EST) of [17] projected onto the helical basis. The solution behaviour of the HST model will be compared to that for the EST model for several data assimilation objectives formulated on the helical triad modes. For classical deterministic models one obtains a system of ordinary differential equations. However, for stochastic dynamics a set of stochastic differential equations (SDEs) is obtained [8, 14, 25].

The goal of the data assimilation procedure in this context is twofold: firstly, it is used to calibrate the uncertainty of the model (the amplitude of the noise). Secondly, once the calibration is complete, the particle filtering methodology can be used to reduce the uncertainty. We want the distribution of the fluctuations to be properly approximated. In the absence of stochasticity, all particles would go in the same direction and the initial spread would rapidly disappear because of the hyperbolic character of the model. In the absence of a reasonable spread, the particle filter methodology will eventually collapse. For this reason, we need to introduce stochasticity into the system that correctly characterises the fluctuation dynamics. In particular, we want to find the type of noise amplitude and the stochastic parameters for which the distribution of the output samples is reasonably uniform.

**Structure of the Paper** In Sect. 2 we introduce the triad models for incompressible flows modelled by the Euler equation in its deterministic and stochastic form. To this end, we introduce the stochastic parametrisation paradigms. Building upon these models for the 3D Euler equation, we then present reduced order triad models derived from the original equations. The derivation follows the classical approach for triad models from the literature, that has been successfully employed in the deterministic case. Our full derivations, complete also for the stochastically parametrised models, can be found in Appendix 2. The Data Assimilation experiments are carried out in Sect. 3. We first briefly explain the standard particle filter methodology, and then in Sect. 3.1 we present the findings of our numerical studies. In particular, Sect. 3.1 presents the results of the main numerical studies in this work. These are:
* The model realisations of the stochastic models which we propose in this work are presented in Sect. 3.1.3 for different realisations of the noise. Further, we numerically confirm here the physical conservation properties of the models that motivated their theoretical design.

* The model statistics for a large number of stochastic realisations are presented in Sect. 3.1.4. Crucially, the evolution of the model statistics in time reveals differences between the two stochastic models on the level of the individual triad energies which go beyond the conservation properties. Moreover, we show here the (in-)stability of either model with respect to large noise parameters.

* In Sect. 3.1.5 we exhibit the results of the data assimilation performed using the transition kernels derived from the theoretical models developed in this paper. We show here that using the stochastic transition kernels associated with our proposed stochastic models improves the particle filtering procedure and produces efficient ensemble evolutions that are well-suited for data assimilation purposes.

In Sect. 4 we describe our conclusions on this topic. We conclude the paper with a number of appendices: in Appendix 1 one can find a list of notations and standard identities, in Appendix 2 we present a detailed derivation of shell models (deterministic and stochastic), in Appendix 3 we introduce some supplementary numerics related to the noise amplitude calibration.

**Code Availability** The code corresponding to the numerical experiments in this paper is archived in [22]. The GitHub repository is located at https://github.com/alobbe/stochastic-triads.

# 2 Reduced Order Models for Incompressible Fluids

## 2.1 Reduced Order Models for the 3D Euler Equation

The 3D Euler equations model incompressible inviscid fluid dynamics. These equations may be written by using the Leray operator $\mathcal{P}$ to project onto the divergence-free part of its operand as

$$ \begin{aligned} \frac{\partial \mathbf{v}}{\partial t} &= \mathcal{P}(\mathbf{v} \times \mathrm{curl} \mathbf{v}) \\ &= \mathcal{P}\left(\frac{\delta E}{\delta \mathbf{v}} \times \frac{\delta C}{\delta \mathbf{v}}\right) \end{aligned} \tag{1} $$

$$ \text{with conserved Energy } E(\mathbf{v}) = \int_{\mathbb{R}^3} \frac{1}{2} \mathcal{P} \mathbf{v} \cdot \mathbf{v} d^3 x $$

$$ \text{and conserved helicity } C(\mathbf{v}) = \frac{1}{2} \int_{\mathbb{R}^3} \mathbf{v} \cdot \mathrm{curl} \mathbf{v} d^3 x $$

where $\delta/\delta\mathbf{v}$ represents variational derivative with respect to the fluid velocity $\mathbf{v}$.
Following [9, 10] we use a Galerkin expansion in orthogonal vector modes that are eigenfunctions of the curl operator. Assume the fluid is contained in a periodic box $D \subset \mathbb{R}^3$ of side length $L > 0$. Then the velocity $\mathbf{v}$ and vorticity $\omega := \nabla \times \mathbf{v}$ may be expanded in circularly polarised or helical modes $\mathbf{h}_{\pm}(\mathbf{k}) \exp(i\mathbf{k} \cdot \mathbf{x})$, with the wave vectors $\mathbf{k} \in \mathcal{K} := (2\pi/L)\mathbb{Z}^3$. The modes shall be orthogonal, i.e.

$$ \int_D \mathbf{h}_{s_p}(\mathbf{p}) \exp(i\mathbf{p} \cdot \mathbf{x}) \cdot [\mathbf{h}_{s_q}(\mathbf{q}) \exp(i\mathbf{q} \cdot \mathbf{x})]^* \mathrm{d}\mathbf{x} = C \delta_{\mathbf{p},\mathbf{q}} \delta_{s_p,s_q}; \quad C > 0 \text{ const.} \quad (2) $$

The complex vector amplitudes $\mathbf{h}_{\pm}(\mathbf{k})$ should satisfy $\mathbf{k} \cdot \mathbf{h}_{\pm}(\mathbf{k}) = 0$ and $i\mathbf{k} \times \mathbf{h}_{\pm}(\mathbf{k}) = \pm |\mathbf{k}|\mathbf{h}_{\pm}(\mathbf{k})$. A convenient choice of basis for the $\mathbf{h}_{\pm}(\mathbf{k})$ is then given by

$$ \mathbf{h}_{\pm}(\mathbf{k}) := \nu \times \kappa \pm i\nu, \quad \text{with } \kappa := \mathbf{k}/k, \quad \nu := \mathbf{k} \times \mathbf{\Gamma} / |\mathbf{k} \times \mathbf{\Gamma}|, \quad \mathbf{\Gamma} := \text{const}, \quad (3) $$

for which $|\mathbf{h}_{\pm}(\mathbf{k})|^2 := \mathbf{h}_{\pm}(\mathbf{k}) \cdot \mathbf{h}_{\pm}(\mathbf{k})^* = 2$ and $\mathbf{h}_{\pm}(\mathbf{k}) \cdot \mathbf{h}_{\mp}(\mathbf{k})^* = 0$.

At this point, one notices the key features of the helical modes $\mathbf{h}_{\pm}(\mathbf{k}) \exp(i\mathbf{k} \cdot \mathbf{x})$ which greatly simplifies analysis of modal expansions of the 3D Euler and related equations, such as 3D Navier-Stokes. Namely, the helical modes $\mathbf{h}_{\pm}(\mathbf{k}) \exp(i\mathbf{k} \cdot \mathbf{x})$, are divergence-free eigenfunctions of the curl operator. Specifically,

$$ \nabla \cdot \mathbf{h}_s(\mathbf{k}) e^{i\mathbf{k} \cdot \mathbf{x}} = i\mathbf{k} \cdot \mathbf{h}_s(\mathbf{k}) e^{i\mathbf{k} \cdot \mathbf{x}} = 0 \quad (4) $$

and

$$ \nabla \times \mathbf{h}_s(\mathbf{k}) e^{i\mathbf{k} \cdot \mathbf{x}} = i\mathbf{k} \times \mathbf{h}_s(\mathbf{k}) e^{i\mathbf{k} \cdot \mathbf{x}} = s|\mathbf{k}|\mathbf{h}_s(\mathbf{k}) e^{i\mathbf{k} \cdot \mathbf{x}}. \quad (5) $$

See [9, 10, 28, 29] for more information about how this Galerkin decomposition into divergence-free eigenfunctions of the curl are used as a standard tool in analysis of 3D solution behaviour of the deterministic Euler fluid equations and Navier-Stokes fluid equations. In particular, the helical mode expansions in Eqs. (6) comprise the source of the popular shell models as finite-dimensional expansions of turbulent fluid dynamics. Thus, this expansion provides a useful framework for studying low-dimension stochastic models of 3D Navier-Stokes turbulence.

In terms of the basis of helical modes $\mathbf{h}_{\pm}(\mathbf{k}) \exp(i\mathbf{k} \cdot \mathbf{x})$, the divergence-free fluid velocity $\mathbf{v}(\mathbf{x}, t)$ and vorticity $\omega(\mathbf{x}, t)$ are expressed in [28, 29] in terms of complex vector amplitudes $\mathbf{u}(\mathbf{k}, t)$, $\mathbf{\omega}(\mathbf{k}, t) \in \mathbb{C}^3$, respectively,

$$ \begin{aligned} \mathbf{v}(\mathbf{x}, t) &:= \sum_{\mathbf{p}} \mathbf{u}(\mathbf{p}, t) e^{i\mathbf{p} \cdot \mathbf{x}} := \sum_{\mathbf{p}} \sum_{s_p=\pm} a_{s_p}(\mathbf{p}, t) \mathbf{h}_{s_p}(\mathbf{p}) e^{i\mathbf{p} \cdot \mathbf{x}}, \\ \omega(\mathbf{x}, t) &:= \sum_{\mathbf{q}} \mathbf{\omega}(\mathbf{q}, t) e^{i\mathbf{q} \cdot \mathbf{x}} := \sum_{\mathbf{q}} \sum_{s_q=\pm} s_q |\mathbf{q}| a_{s_q}(\mathbf{q}, t) \mathbf{h}_{s_q}(\mathbf{q}) e^{i\mathbf{q} \cdot \mathbf{x}}. \end{aligned} \quad (6) $$
Here, the choice

$$ \mathbf{u}(\mathbf{k}, t) := a_+(\mathbf{k}, t)\mathbf{h}_+(\mathbf{k}) + a_-(\mathbf{k}, t)\mathbf{h}_-(\mathbf{k}) = \sum_{s_k=\pm} a_{s_k}(\mathbf{k}, t)\mathbf{h}_{s_k}(\mathbf{k}), \quad (7) $$

with $a_s^*(\mathbf{k}) = a_s(-\mathbf{k})$ [29], was made, so that (5) implies

$$ \bar{\omega}(\mathbf{k}, t) := |\mathbf{k}| (a_+(\mathbf{k}, t)\mathbf{h}_+(\mathbf{k}) - a_-(\mathbf{k}, t)\mathbf{h}_-(\mathbf{k})) = |\mathbf{k}| \sum_{s_k=\pm} s_k a_{s_k}(\mathbf{k}, t)\mathbf{h}_{s_k}(\mathbf{k}). \quad (8) $$

The conservation laws for the Euler fluid kinetic energy and helicity—expressed as integrals over the spatially periodic box $D$—can be evaluated in Fourier space via Parseval’s theorem, as follows,

$$ \begin{aligned} \frac{1}{2} \int_D |\mathbf{v}(\mathbf{x}, t)|^2 d^3x &= \frac{1}{2} \sum_{\mathbf{k}} \mathbf{u}(\mathbf{k}, t) \cdot \mathbf{u}^*(\mathbf{k}, t) = \sum_{\mathbf{k}} \sum_{s_k=\pm} a_{s_k}(\mathbf{k}, t) \cdot a_{s_k}^*(\mathbf{k}, t), \\ \int_D \mathbf{v}(\mathbf{x}, t) \cdot \mathrm{curl} \mathbf{v}(\mathbf{x}, t) d^3x &= \sum_{\mathbf{k}} \mathbf{u}(\mathbf{k}, t) \cdot \bar{\omega}^*(\mathbf{k}, t) \\ &= \sum_{\mathbf{k}} \sum_{s_k=\pm} k s_k a_{s_k}(\mathbf{k}, t) a_{s_k}^*(\mathbf{k}, t) \mathbf{h}_{s_k}(\mathbf{k}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \\ &= 2 \sum_{\mathbf{k}} \sum_{s_k=\pm} k s_k a_{s_k}(\mathbf{k}, t) a_{s_k}^*(\mathbf{k}, t). \end{aligned} \quad (9) $$

Expanding the terms of the Euler equation in curl form (57), we obtain the Euler equations for the coefficients $a_{s_k}(\mathbf{k}, t)$. For all $\mathbf{k} \in \mathcal{K}$, $s_k \in \{+, -\}$ we have

$$ \partial_t a_{s_k}(\mathbf{k}, t) = -\frac{1}{4} \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} (s_p |\mathbf{p}| - s_q |\mathbf{q}|) a_{s_p}^*(\mathbf{p}, t) a_{s_q}^*(\mathbf{q}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}). \quad (10) $$

For the explicit derivation of Eq. (10) see section “Deterministic Euler” in Appendix 2.

The elementary interactions in Fourier space take place between triads of wave vectors such that $\mathbf{k} + \mathbf{p} + \mathbf{q} = 0$, as exemplified in Eq. (10) above. There are two degrees of freedom per wave vector, $(a_+, a_-)$, so eight different types of interaction are allowed according to the value of the triplet $(s_k, s_p, s_q) = (\pm 1, \pm 1, \pm 1)$. Consider a fixed triple of wave vectors $\mathbf{k}, \mathbf{p}, \mathbf{q} \in \mathcal{K}$ such that $\mathbf{k} + \mathbf{p} + \mathbf{q} = 0$ and a fixed triple $s_k, s_p, s_q \in \{+, -\}$. This gives rise to three coefficients $a_{s_k}(\mathbf{k}, t), a_{s_p}(\mathbf{p}, t), a_{s_q}(\mathbf{q}, t)$, which we compactly summarise into the complex vector
$$ \mathbf{a} = (a_{s_k}, a_{s_p}, a_{s_q}) \in \mathbb{C}^3. $$

The dynamics of **a** is determined by the three equations obtained from (10)

$$ \frac{da_{s_k}}{dt} = g(s_p|\mathbf{p}| - s_q|\mathbf{q}|)a_{s_p}^* a_{s_q}^* + R, \quad (11) $$

$$ \frac{da_{s_p}}{dt} = g(s_q|\mathbf{q}| - s_k|\mathbf{k}|)a_{s_q}^* a_{s_k}^* + R, \quad (12) $$

$$ \frac{da_{s_q}}{dt} = g(s_k|\mathbf{k}| - s_p|\mathbf{p}|)a_{s_k}^* a_{s_p}^* + R, \quad (13) $$

We pick out a summand and cycle through $k, p, q$

$$ \frac{da_{s_k}}{dt} = g(s_p|\mathbf{p}| - s_q|\mathbf{q}|)a_{s_p}^* a_{s_q}^*, \quad (14) $$

$$ \frac{da_{s_p}}{dt} = g(s_q|\mathbf{q}| - s_k|\mathbf{k}|)a_{s_q}^* a_{s_k}^*, \quad (15) $$

$$ \frac{da_{s_q}}{dt} = g(s_k|\mathbf{k}| - s_p|\mathbf{p}|)a_{s_k}^* a_{s_p}^*, \quad (16) $$

with the constant complex scalar

$$ g := -\frac{1}{4} \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}). \quad (17) $$

The equations corresponding to the single triad interaction of type $(s_k, s_p, s_q)$ with $\mathbf{k} + \mathbf{p} + \mathbf{q} = 0$ thus have the complex vector form also derived in [28],

$$ \frac{da}{dt} = g\mathbf{a}^* \times \mathbf{D}\mathbf{a}^* = g(\mathbf{a} \times \mathbf{D}\mathbf{a})^*, \quad (18) $$

with the constant diagonal matrix

$$ \mathbf{D} := \mathrm{diag}(s_k|\mathbf{k}|, s_p|\mathbf{p}|, s_q|\mathbf{q}|). \quad (19) $$

The form of the factor $g$ defined in (17) above can be calculated from (3) to show that it depends on the shape and the orientation of the wave-vector triad, but not on its scale; since the real and imaginary parts of the complex helical vector amplitudes $\mathbf{h}_{s_k}$ are unit vectors. Moreover, $\mathbf{D}\mathbf{a}$ can be seen to represent the $(s_k, s_p, s_q)$ components of the vorticity vector amplitude $\boldsymbol{\omega}$ through Eq. (8) above. Two conservation laws for real-valued triad energy and helicity follow immediately from Eq. (18), as

$$ \frac{d}{dt}(\mathbf{a} \cdot \mathbf{a}^*) = 0 \quad \text{and} \quad \frac{d}{dt}(\mathbf{a} \cdot \mathbf{D}\mathbf{a}^*) = 0. \quad (20) $$
The dynamical system in Eq. (18) is similar to rigid body dynamics, but replaced by complex angular momentum (a), complex angular velocity (Da) and real moment of inertia I = D⁻<fnref n="1" /> relating the two complex quantities.<fnref n="1" />

## 2.2 Stochastic Parametrizations for the 3D Euler Equation

In the following we introduce a reduced order model for the stochastic parametrizations introduced through the Stochastic Advection by Lie Transport paradigm as well as the Location Uncertainty paradigm. We explain below the rationale of these parametrizations:

### 2.2.1 Modelling Under the Stochastic Advection by Lie Transport Principle

The SALT equations were derived in [16] using a Stratonovich stochastic version of Hamilton's variational principle [18] in combination with Kraichnan's scalar turbulence model based on Stratonovich stochastic Lagrangian paths [20]. The application of Hamilton's principle with an imposed stochastic Lie transport constraint implied an Euler-Poincaré equation for the fluid motion [18]. The 3D SALT Euler equations for divergence-free fluid velocity $v(x, t)$ are given by

$$ 
\begin{gathered}
d\mathbf{v} + (d\mathbf{x}_t \cdot \nabla)\mathbf{v} + v_j \nabla \mathbf{x}_t^j = -\nabla dp, \\
\text{with } d\mathbf{x}_t = \mathbf{v}dt + \sum_i \xi_i(x) \circ dW_t^i.
\end{gathered}
\quad (21) $$

As discussed in [16], this motion equation yields a Kelvin-Noether circulation theorem for the stochastic system

$$ d \oint_{c(x_t)} \mathbf{v} \cdot d\mathbf{x} = - \oint_{c(x_t)} \nabla dp \cdot d\mathbf{x} = 0. \quad (22) $$

This stochastic Kelvin circulation theorem has the same form as that for the deterministic system, except that each line element of the material loop in Kelvin's theorem follows the Stratonovich stochastic Lagrangian path, $x_t$.

The real vectors $\xi_i$ comprise the time-independent noise amplitudes which are to be determined from data assimilation. The $W^i$ are independent (uni-dimensional) standard Brownian motions and $\circ$ denotes stochastic integration in the Stratonovich sense.<fnref n="2" /> The curl form of the SALT Euler motion equation in (21) is obtained

<footnote n="1">¹ Rigid body dynamics with complex angular momentum has also been discussed previously in [3].</footnote>

<footnote n="2">² An exposition of Brownian motion, stochastic calculus and the Stratonovich integral is to be found, for example, in the monograph [19].</footnote>
from (56) and given by

$$d\mathbf{v} - d\mathbf{x}_t \times \mathrm{curl}\mathbf{v} + \nabla(\mathbf{v} \cdot d\mathbf{x}_t) = -\nabla dp. \quad (23)$$

The motion equation (23) and its curl yielding the SALT vorticity equation implies a formula for the evolution of the helicity of the flow, $\Lambda$, defined as

$$\Lambda := \int_D \mathbf{v} \cdot \mathrm{curl}\mathbf{v} \, d^3x. \quad (24)$$

Upon applying the divergence theorem, one finds

$$d\Lambda = - \int_{\partial D} \hat{\mathbf{n}} \cdot \left( (\mathbf{v} \cdot \mathrm{curl}\mathbf{v}) d\mathbf{x}_t + \mathrm{curl}\mathbf{v} \, dp \right) dS. \quad (25)$$

Thus, in a periodic 3D domain, or in an infinite 3D domain with asymptotically vanishing boundary conditions, the SALT motion equation in (21) or (23) preserves the helicity, $\Lambda$, defined in (24). However, a glance at the SALT motion equation in (23) informs us that it will not preserve the kinetic energy, since even with the usual fluid boundary conditions $\mathrm{div}\mathbf{v} = 0$ implies

$$\frac{1}{2} d\|\mathbf{v}\|_{L^2}^2 := \int_D \mathbf{v} \cdot d\mathbf{x}_t \times \mathrm{curl}\mathbf{v} \, d^3x \neq 0. \quad (26)$$

### 2.2.2 Modeling Under the Location Uncertainty Principle

The Location Uncertainty principle consists in decomposing the flow trajectory $\mathbf{x}: \Omega \times \mathbb{R}^+ \to \Omega$ over a bounded domain, $\Omega \subset \mathbb{R}^3$

$$d\mathbf{x}_t = \mathbf{v}(\mathbf{x}_t, t) \, dt + \boldsymbol{\sigma}(\mathbf{x}_t, t) \, d\mathbf{W}_t \quad (27)$$

in terms of $\mathbf{v}(\mathbf{x}_t, t)$, a smooth-in-time component of the (Lagrangian) velocity and a noise $\boldsymbol{\sigma}(\mathbf{x}_t, t) d\mathbf{W}_t$, which has here to be understood in the Itô sense and that accounts for the unresolved processes. The Wiener process, $\mathbf{W}_t$ is a $H$-valued (cylindrical) Brownian motion, where $H$ is the Hilbert space of square integrable functions. The noise is then properly defined as the application of an Hilbert-Schmidt symmetric integral kernel $\boldsymbol{\sigma}_t f(\mathbf{x}) = \int_S \check{\boldsymbol{\sigma}}(\mathbf{x}, \mathbf{y}, t) f(\mathbf{y}) \, d\mathbf{y}$ to the $H$-valued cylindrical Wiener process $\mathbf{W}$ as

$$(\boldsymbol{\sigma}_t d\mathbf{W}_t)^i(\mathbf{x}) = \int_S \check{\boldsymbol{\sigma}}_{ik}(\mathbf{x}, \mathbf{y}, t) dW_t^k(\mathbf{y}) \, d\mathbf{y}, \quad (28)$$
The role of the integrable kernel $\check{\sigma}$ is to impose a spatial correlation on the small-scale component. It leads to the covariance tensor $Q$

$$Q_{ij}(\mathbf{x}, \mathbf{y}, t, s) = \mathrm{E} \left[ \left( \boldsymbol{\sigma}_t \mathrm{d}\mathbf{W}_t(\mathbf{x}) \right)^i \left( \boldsymbol{\sigma}_t \mathrm{d}\mathbf{W}_s(\mathbf{y}) \right)^j \right] \\ = \delta(t-s) \mathrm{d}t \int_S \check{\sigma}_{ik}(\mathbf{x}, \mathbf{z}, t) \check{\sigma}_{kj}(\mathbf{z}, \mathbf{y}, s) \mathrm{d}\mathbf{z},$$

of the centered Gaussian process $\boldsymbol{\sigma}_t \mathrm{d}\mathbf{W}_t \sim \mathcal{N}(0, \mathbf{Q}dt)$. The diagonal components of the covariance tensor per unit of time, referred to as the variance tensor, $\mathbf{a}$, is a positive definite matrix defined as $\mathbf{a}(\mathbf{x}, t)\delta(t-t')\mathrm{d}t = \mathbf{Q}(\mathbf{x}, \mathbf{x}, t, t')$, that quantifies the strength of the noise and has the dimension of a viscosity in $\mathrm{m}^2\mathrm{s}^{-1}$. The operator $\mathbf{Q}$ being compact auto-adjoint positive definite operator on $H$, it admits eigenfunctions $\xi_n(\cdot, t)$ with (strictly) positive eigenvalues $\lambda_n(t)$ satisfying $\sum_{n \in \mathbb{N}} \lambda_n(t) < +\infty$. As a consequence, the noise and the variance tensor $\mathbf{a}$ can be expressed through the spectral representation

$$\boldsymbol{\sigma}_t \mathrm{d}\mathbf{W}_t(\mathbf{x}) = \sum_{n \in \mathbb{N}} \lambda_n^{1/2}(t) \xi_n(\mathbf{x}, t) \mathrm{d}\beta_n \quad (29)$$

$$\mathbf{a}(\mathbf{x}, t) = \sum_{n \in \mathbb{N}} \lambda_n(t) \xi_n(\mathbf{x}, t) \xi_n^{\dagger}(\mathbf{x}, t). \quad (30)$$

The rate of change of a volume $V_t$ of the scalar $q$ is given by the stochastic Reynolds transport theorem, introduced in [23]

$$\mathrm{d} \int_{V_t} q(\mathbf{x}, t) \mathrm{d}\mathbf{x} = \int_{V_t} \{ \mathrm{D}_t q + q \nabla \cdot [\mathbf{v}^* \mathrm{d}t + \boldsymbol{\sigma}_t \mathrm{d}\mathbf{W}_t] \} (\mathbf{x}, t) \mathrm{d}\mathbf{x}, \quad (31)$$

with the transport operator

$$\mathrm{D}_t q = \mathrm{d}_t q + [\mathbf{v}^* \mathrm{d}t + \boldsymbol{\sigma}_t \mathrm{d}\mathbf{W}_t] \cdot \nabla q - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla q) \mathrm{d}t. \quad (32)$$

In this formula, the first component of the right-hand side is the increment in time at a fixed location of the process $q$, that is $\mathrm{d}_t q = q(\mathbf{x}_t, t+\mathrm{d}t) - q(\mathbf{x}_t, t)$, playing the role of a derivative in time for a non differentiable process. The effective velocity $\mathbf{v}^*$ is defined as

$$\mathbf{v}^* = \mathbf{v} - \frac{1}{2} \nabla \cdot \mathbf{a} + \boldsymbol{\sigma}_t^* (\nabla \cdot \boldsymbol{\sigma}_t), \quad (33)$$

where the velocity component $\mathbf{v}_s = \frac{1}{2} \nabla \cdot \mathbf{a}$ results from the noise inhomogeneities. For incompressible homogeneous noise as considered in this work $\mathbf{v}^* = \mathbf{v}$. Besides, the diffusion term exactly balances the noise brought by the noise. With Stratonovich convention and a homogeneous noise the transport operator takes a
simplified form similar to the material derivative:

$$D_t q = d_t q + (\mathbf{v} dt + \boldsymbol{\sigma}_t \circ d\mathbf{W}_t) \cdot \nabla q. \quad (34)$$

For a divergence free homogeneous noise, the Euler equation, in the LU framework, can then be defined as:

$$d_t \mathbf{v} + (\mathbf{v} dt + \boldsymbol{\sigma}_t \circ d\mathbf{W}_t) \cdot \nabla \mathbf{v} = -\nabla d_p t, \quad \nabla \cdot \mathbf{v} = 0, \quad (35)$$

where $d_p t$ denotes the pressure composed of a finite variation term and a martingale pressure term. With the Leray projection, $\mathbb{P}$, this pressure term can be removed and we obtain, the inertial form of the Euler equation:

$$d_t \mathbf{v} + \mathbb{P}((d\mathbf{x}_t \cdot \nabla)\mathbf{v}) = 0, \quad \nabla \cdot \mathbf{v} = 0. \quad (36)$$

## 2.3 Triad Model Comparison

The reduced order model for the full-scale 3D SALT Euler and 3D LU Euler for a single triad interaction equation is obtained by projecting the continuous stochastic Euler models onto the helical modes, in the same fashion as for the deterministic equation (18). Therefore, we introduce an additional Stratonovich stochastic term into the transport velocity in (7) as

$$d\mathbf{x}_t(\mathbf{k}, t) := (a_+(\mathbf{k}, t)\mathbf{h}_+ + a_-(\mathbf{k}, t)\mathbf{h}_-) dt + \sum_i (b_+^i(\mathbf{k})\mathbf{h}_+ + b_-^i(\mathbf{k})\mathbf{h}_-) \circ dW_t^i, \quad (37)$$

where the $\mathbf{k}$-dependent complex vector $\mathbf{b}(\mathbf{k}) := (b_{s_k}, b_{s_p}, b_{s_q})^T \in \mathbb{C}^3$ represents the time-independent noise amplitude which is to be determined from data assimilation, similar to the continuous stochastic models (21). Enumerating the equation for a single triad then yields, after rearranging using exchange symmetry in $(\mathbf{k}, \mathbf{p}, \mathbf{q})$, the matrix equation

$$d \begin{bmatrix} a_{s_k} \\ a_{s_p} \\ a_{s_q} \end{bmatrix} = g \begin{bmatrix} 0 & -q_{s_q} a_{s_q} & p_s a_{s_p} \\ q_{s_q} a_{s_q} & 0 & -k_{s_k} a_{s_k} \\ -p_s a_{s_p} & k_{s_k} a_{s_k} & 0 \end{bmatrix}^* \begin{bmatrix} a_{s_k} \, dt + b_{s_k} \circ dW_t \\ a_{s_p} \, dt + b_{s_p} \circ dW_t \\ a_{s_q} \, dt + b_{s_q} \circ dW_t \end{bmatrix}^*. \quad (38)$$

Upon applying the previous steps for the deterministic case to the stochastic velocity in (37), the single triad interaction dynamics for the SALT case would emerge as, cf. Eq. (18),

$$d\mathbf{a} = g(\mathbf{a}(\mathbf{k}, t)dt + \mathbf{b}(\mathbf{k}) \circ dW_t)^* \times \mathbb{D}\mathbf{a}^*. \quad (39)$$
The details of this computation can be found in section “LU Euler” in Appendix 2 for the 3D LU Euler model and in section “SALT Euler” in Appendix 2 for the 3D SALT Euler model.

Remarkably, the HST equation for triad interaction (39) still preserves the triad helicity $\mathbf{a} \cdot \mathbf{D}\mathbf{a}^*$. Hence we name this model the *helicity preserving stochastic triad* (HST) model. Note that in both equations we use as single source of noise (One Brownian motion drives the entire system).

It is readily checked that the HST triad evolution (39) preserves the helicity. Let's have a look at the diffusion coefficients.

$$ \mathbf{a}^* \times \mathbf{D}\mathbf{b} = \begin{bmatrix} a_p^* s_q q b_q - a_q^* s_p p b_p \\ a_q^* s_k k b_k - a_k^* s_q q b_q \\ a_k^* s_p p b_p - a_p^* s_k k b_k \end{bmatrix}, \quad \mathbf{b} \times \mathbf{D}\mathbf{a}^* = \begin{bmatrix} b_p s_q q a_q^* - b_q s_p p a_p^* \\ b_q s_k k a_k^* - b_k s_q q a_q^* \\ b_k s_p p a_p^* - b_p s_k k a_k^* \end{bmatrix}. \quad (40) $$

Taking the difference

$$ \mathbf{a}^* \times \mathbf{D}\mathbf{b} - \mathbf{b} \times \mathbf{D}\mathbf{a}^* = \begin{bmatrix} a_p^*(s_q q + s_p p)b_q - a_q^*(s_p p + s_q q)b_p \\ a_q^*(s_k k + s_q q)b_k - a_k^*(s_q q + s_k k)b_q \\ a_k^*(s_p p + s_k k)b_p - a_p^*(s_k k + s_p p)b_k \end{bmatrix}. \quad (41) $$

Writing $\rho := \mathrm{Tr} \mathbf{D}$ we get

$$ \mathbf{a}^* \times \mathbf{D}\mathbf{b} - \mathbf{b} \times \mathbf{D}\mathbf{a}^* = \begin{bmatrix} a_p^*(\rho - s_k k)b_q - a_q^*(\rho - s_k k)b_p \\ a_q^*(\rho - s_p p)b_k - a_k^*(\rho - s_p p)b_q \\ a_k^*(\rho - s_q q)b_p - a_p^*(\rho - s_q q)b_k \end{bmatrix}. \quad (42) $$

So that the difference term becomes

$$ \mathbf{a}^* \times \mathbf{D}\mathbf{b} - \mathbf{b} \times \mathbf{D}\mathbf{a}^* = (\rho \mathrm{Id} - \mathbf{D})(\mathbf{a} \times \mathbf{b})^*. \quad (43) $$

Since the projections of the LU and SALT models onto a single triad are indistinguishable, we introduce a different model that conserves energy on a single triad to enable a comparison between energy conserving and helicity conserving models.

**Energy-Preserving Stochastic Triad (EST) Model** We introduce below a modified version of the HST triad equation (39) that introduces stochasticity into the vorticity instead of into the transport velocity and thereby conserves the energy. This is inspired by the full-scale model introduced in [17]. The reduced model is as follows

$$ \mathbf{d}\mathbf{a} = -g\mathbf{a}^* \times \mathbf{D}(\mathbf{a}(\mathbf{k}, t) \mathrm{d}t + \mathbf{b}(\mathbf{k}) \circ \mathrm{d}W_t)^*. \quad (44) $$
We call this model the energy preserving stochastic triad (EST). Written in matrix form Eq. (44) becomes

$$ d \begin{bmatrix} a_{s_k} \\ a_{s_p} \\ a_{s_q} \end{bmatrix} = g \begin{bmatrix} 0 & -a_{s_q} & a_{s_p} \\ a_{s_q} & 0 & -a_{s_k} \\ -a_{s_p} & a_{s_k} & 0 \end{bmatrix}^* \begin{bmatrix} k a_k (a_{s_k} dt + b_{s_k}(\mathbf{k}) \circ dW_t) \\ p s_p (a_{s_p} dt + b_{s_p}(\mathbf{p}) \circ dW_t) \\ q s_q (a_{s_q} dt + b_{s_q}(\mathbf{q}) \circ dW_t) \end{bmatrix}^* . \quad (45) $$

The exchange symmetry between the two models HST and EST in the placement of the noise in Eqs. (39) and (44) is apparent already in the exchange symmetry between velocity and vorticity in Euler's fluid equations (1).

**Deviation from the Conservation Laws** We can write the equations for the deviation from the conservation laws, which is present in both models. The SALT model deviates from the energy conservation by

$$ d_t E_{\text{HST}} = g \mathbf{b} \cdot (\mathbf{D}\mathbf{a}^* \times \mathbf{a}^*) \circ dW_t \quad (46) $$

whereas the LU model deviates from the helicity conservation by

$$ d_t H_{\text{EST}} = g \mathbf{D}\mathbf{b} \cdot (\mathbf{D}\mathbf{a}^* \times \mathbf{a}^*) \circ dW_t. \quad (47) $$

This is seen, since, to get the energy we dot the HST equation with $\mathbf{a}^*$ and to get helicity we dot the EST equation with $\mathbf{D}\mathbf{a}^*$ and use the standard identities

$$ \mathbf{a}^* \cdot (\mathbf{b} \times \mathbf{D}\mathbf{a}^*) = \mathbf{b} \cdot (\mathbf{D}\mathbf{a}^* \times \mathbf{a}^*) \quad (48) $$

$$ \mathbf{D}\mathbf{a}^* \cdot (\mathbf{a}^* \times \mathbf{D}\mathbf{b}) = \mathbf{D}\mathbf{b} \cdot (\mathbf{D}\mathbf{a}^* \times \mathbf{a}^*). \quad (49) $$

Therefore, $\mathbf{b}$ respects the right scaling and no further scale adjustments between the SALT and LU noise scaling need to be performed in order to compare the models.

## 3 Data Assimilation Comparison

In this section, we perform a comparative study of the two reduced order models (HST and EST) introduced above by using data assimilation tools. The particular methodology that we make use of is that of particle filters. We will first briefly explain the particle filtering methodology in a generic framework:

Let $X$ and $Z$ be two processes defined on a given probability space $(\Omega, \mathcal{F}, \mathbb{P})$. The process $X$ is usually called the *signal process* or the *truth* and $Z$ is the *observation process*. In this paper, $X$ is the pathwise solution of a (deterministic) shell model. The pair of processes $(X, Z)$ forms the basis of the nonlinear filtering problem which consists in finding the best approximation of the posterior
distribution of the signal $X_t$ given the observations $Z_1, Z_2, \dots, Z_{t_n}$.<fnref n="3" /> The posterior distribution of the signal at time $t$ is denoted by $\pi_t$. We let $d_X$ be the dimension of the state space and $d_Z$ be the dimension of the observation space. This mixed continuous-discrete time framework can be embedded into a fully discrete framework, whereby one is interested in computing the conditional probability law of the signal at the time corresponding to the observation time. in other words one wants to compute the conditional distribution $\pi_{t_n}$ of $X(t_n)$ given the data $Z(t_1), Z(t_2), \dots, Z(t_n)$. The process $X$ is assumed to be a Markov process, and we will denote by $\mathcal{K}_n$ its transition kernel, that is

$$ \mathcal{K}_n : \mathbb{R}^{d_X} \times \mathcal{B}(\mathbb{R}^{d_X}) \to [0, 1], \quad \mathcal{K}_n(x, B) = \mathbb{P}(X_{t_n} \in B | X_{t_{n-1}} = x) \quad (50) $$

for any Borel measurable set $B \in \mathcal{B}(\mathbb{R}^{d_X})$ and $x \in \mathbb{R}^{d_X}$. The process $Z$ models noisy measurements of the truth, using the so-called observation operator $\mathcal{H}: \mathbb{R}^{d_X} \to \mathbb{R}^{d_Z}$:

$$ Z_n = \mathcal{H}(X_{t_n}) + V_n \quad (51) $$

where $\{V_n\}_{n \ge 0}$ are independent identically distributed random variables that represent the measurement noise and $\mathcal{H}$ is a Borel-measurable function. In this paper we will assume that $\{V_n\}_{n \ge 0}$ have standard normal distributions, but the same methodology can be applied to more general distributions. Observations are incorporated into the system at *assimilation times*. The following recursion formula holds (see [2])

$$ \pi_n = g_n \star \pi_{n-1} \mathcal{K}_n \quad (52) $$

where by ‘★’ we denoted the projective product (see e.g. Definition 10.4 in [2]).

In the following, we compare approximations of the posterior distribution of the signal using *particle filters*. These are sequential Monte Carlo methods which generate approximations of the posterior distribution $\pi_t$ using sets of *particles*. That is, they generate approximations that are (random) measures of the form

$$ \pi_n \approx \sum_{\ell} w_n^{\ell} \delta(x_n^{\ell}), $$

where $\delta$ is the Dirac delta distribution, $w_t^1, w_t^2, \dots$ are the *weights* of the particles and $x_t^1, x_t^2, \dots$ are their corresponding positions. Particle filters are used to make inferences about the signal process by using Bayes' theorem, the time-evolution induced by the signal $X$, and the observation process $Z$.

<footnote n="3">³ For a mathematical introduction on the subject, see e.g. [2].</footnote>
In a standard particle filter, the particles evolve between assimilation times according to the law of the signal. As we explain below, at each assimilation time the observation is incorporated into the system through the likelihood function:

$$g_t^{zt}: \mathbb{R}^{dx} \to \mathbb{R}_+, g_t^{zt}(x) = g_t(z_t - \mathcal{H}(x)) \quad \text{such that} \quad \mathbb{P}(Z_t \in dz_t | X_t = x) = g_t^{zt}(x) dz_t \quad (53)$$

and all particles are weighted depending on the likelihood of their corresponding position, given the observation. More precisely, the particle $\ell$ is given the weight $w_n^\ell = g_n^{Z^n}(x_\ell)$. Heuristically, the particle weight measures how close the particle trajectory is to the signal trajectory. A selection procedure is then applied to the set of weighted particles. Particles with higher conditional likelihood (guided by the observation) have higher weights and will be multiplied, while those which have small likelihoods will be eliminated. For the basic particle filter, this is done by sampling with replacement from the population of particles, with corresponding probabilities proportional to their weights.

A Monte Carlo implementation of the transition kernel of the signal may not always yield good approximations. In many situations one replaces the original transition kernel with likelihood informed importance proposals, leading to much better approximations. One situation when this is necessary is when the original process is actually deterministic (aside for the initial position which is assumed to be random). This is the case in our paper.

To overcome the collapse of the particle filter when using deterministic transition kernels, one can use a Markov Chain Monte Carlo procedure that leaves the deterministic dynamics invariant. This procedure can be costly and might not always introduce enough spread into the sample. In this paper, we propose a different approach, which we illustrate numerically in Sect. 3.1.2 below. In particular, we propose two different transition kernels based on the physical conservation properties:

* The transition kernel associated with the HST model equation (38). We will denote this transition kernel by $K^{HST}$. As we have explained above, this transition kernel preserves the helicity of the system.

* The transition kernel associated with the triad model equation (45) we will denote this transition kernel by $K^{EST}$. As we have explained above, this transition kernel preserves the energy of the system.

# 3.1 Numerical Studies

## 3.1.1 Numerical Implementation

The models are discretised using the stochastic SSPRK3 scheme which is documented, for example, in [11]. In our specific case, for instance, the HST model is discretised as
$$ \mathbf{q}_1^n = \mathbf{a}_n + g(\mathbf{a}_n^* \times \mathbb{D}\mathbf{a}_n^*) \Delta t + g(\mathbf{b} \times \mathbb{D}\mathbf{a}_n^*) \Delta W $$

$$ \mathbf{q}_2^n = (3/4)\mathbf{a}_n + (1/4)(\mathbf{q}_1^n + g((\mathbf{q}_1^n)^*) \times \mathbb{D}(\mathbf{q}_1^n)^*) \Delta t + g(\mathbf{b} \times \mathbb{D}(\mathbf{q}_1^n)^*) \Delta W $$

$$ \mathbf{a}_{n+1} = (1/3)\mathbf{a}_n + (2/3)(\mathbf{q}_2^n + g((\mathbf{q}_2^n)^*) \times \mathbb{D}(\mathbf{q}_2^n)^*) \Delta t + g(\mathbf{b} \times \mathbb{D}(\mathbf{q}_2^n)^*) \Delta W $$

where $\Delta t$ denotes the timestep and $\Delta W$ the increment of the driving Brownian motion. Further, $\mathbf{a}_n$ is the approximate complex vector amplitude at time $t = n\Delta t$. The EST model is discretised completely analogously. For the numerical simulations we chose the following triad throughout. We set

$$ \mathbf{k} = [1, 0, 0], \quad \mathbf{p} = [0, -1, 1], \quad \mathbf{q} = [-1, 1, -1], \qquad (54) $$

with parities $s_k = 1, s_p = -1, s_q = -1$ and the initial value $\mathbf{a}_0 = \frac{1}{\sqrt{3}}[1, 1, 1]$. We set the parameter $\Gamma = [1, 1, 1]$ and used a time stepsize of $\Delta t = 0.0005$.

### 3.1.2 Data Assimilation for the Deterministic Model

We illustrate the failure of the particle filter with deterministic transition kernel in Fig. 1. In this case, the particle filtering is performed for an ensemble of $n = 25$ particles evolving according to the deterministic triad dynamics. The initial ensemble is spread around the initial value $\mathbf{a}_0$ of the signal according to a Gaussian distribution with standard deviation $1/\sqrt{600}$ and, in particular, does not contain the true initial point. Data assimilation is performed every 10 time units and the observations are taken from the modal energies of the truth with an observation error $\eta \sim \mathcal{N}(0, C)$ with covariance $C = \text{diag}(0.005^2, 0.05^2, 0.05^2)$. We observe that both the bias and the RMSE keep increasing with time to values much larger than the observation error. Moreover, the number of distinct particles decreases rapidly: after 30 steps, a single particle remains that is not the true particle since the true particle was not part of the initial cloud. The particle filter does not work.

### 3.1.3 Reduced Order Model Realisations

The deterministic model in Fig. 1a exhibits continually oscillating triad amplitudes. Plotted are the modal energies. Writing $\mathbf{a} = [a_k, a_p, a_q]$ we call the real value $a_k a_k^*$ the energy of mode $k$. Similarly for the two other modes.

We simulate the model realisations for different noise scenarios. We simulate the effect of noise in each single mode. Let the noise amplitude vector be $\mathbf{b} = [b_k, b_p, b_q]$. Then we simulate the two models for $\mathbf{b} = [b_k = 0.1, b_p = 0, b_q = 0], \mathbf{b} = [b_k = 0, b_p = 0.1, b_q = 0], \text{ and } \mathbf{b} = [b_k = 0, b_p = 0, b_q = 0.1]$. The trajectories of the modal energies for $n = 20$ realisations of the driving noise for each scenario are shown in Fig. 2. We also simulate the case of full noise for
<figure><img src="image_46.png" /><figcaption>**Fig. 1** Deterministic Triad Model. The horizontal axis shows time. (a) Evolution of the modal energies (colored lines) as well as the total energy (dashed line) and helicity (dash-dotted line). The deterministic model exhibits continually oscillating modal energies. The simulation also confirms the conservation of energy and helicity. (b–e) Data assimilation for the deterministic model using particle filter. (b) Evolution of the energy of mode **p** of the signal (grey line) and evolution of the energy of mode **p** for the particle ensemble (blue lines). Noisy observations (black stars) are made and assimilated every 10 time units. (c) The number of unique particle positions in the filtering ensemble. (d) The bias of the particle ensemble wrt. the observations. (e) The RMSE of the particle ensemble wrt. the observations</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>**Fig. 2** Model realisations for both stochastic triad models. Plotted are the modal energies (colored lines), total energy (black line) and helicity (grey line). The respective thick lines are the ensemble means, and the thin lines represent the different stochastic realisations. (a + e) The noise coefficient **b** = [0.1, 0, 0]. (b + f) The noise coefficient **b** = [0, 0.1, 0]. (c + g) The noise coefficient **b** = [0, 0, 0.1]. (d + h) The noise coefficient **b** = [0.1, 0.05, 0.01]</figcaption></figure>
<figure><img src="image_48.png" /><figcaption>**Fig. 3** Evolution of statistical moments of the modal energies for both stochastic models in the full noise case. The statistics are computed pointwise in time from an ensemble of 1000 realisations up to a final time of 150. (a + e) Ensemble mean. (b + f) Ensemble standard deviation. (c + g) Ensemble skew. (d + h) Ensemble kurtosis</figcaption></figure>

the noise amplitude vector **b** = [b_k = 0.1, b_p = 0.05, b_q = 0.01] which was calibrated to the data assimilation objective using the procedure explained in section “Calibration of the Noise Amplitude” in Appendix 3. The ensemble of *n* = 20 realisations of the driving noise in the full noise case is depicted in Fig. 2d and h. In all cases, it can be observed that the mean energy amplitudes of the modes are dampened in both stochastic models. Furthermore, we can experimentally verify the conservation of triad energy for the EST model and the conservation of triad helicity for the HST model.

### 3.1.4 Model Statistics

Figure 3 shows various statistics of the generated ensembles of *n* = 1000 particles for the HST and EST triad models in the full noise case introduced above. We plot the ensemble mean, standard deviation, skew, and kurtosis.

The effect of large noise coefficients is exemplified in Fig. 4. We observe that the HST model explodes whereas the EST model is more tolerant to large noise coefficients, and even in the extreme case, does not become unstable in the mean. The ensemble means are computed from *n* = 500 realisations, using the noise coefficient **b** = [0.0, 1.0, 0.0]. Moreover, to stress the EST model, we also ran the same experiment with a noise coefficient of **b** = [0.0, 10.0, 0.0] for the EST model alone.

The mean ensemble for a large number of particles, *n* = 20,000, is shown in Fig. 5. We can observe that, compared to Fig. 3a and e the oscillations after time 40 are reduced for a very large number of particles. Hence, we believe that the system stabilizes in the mean to stationary modal energies as the limiting effect of the noise.
<figure><img src="image_49.png" /><figcaption>Fig. 4 The effect of large noise coefficients on the mean. Evolution of the mean modal energies (colored lines), mean total energy (black), and mean helicity (grey) for the HST (a) and EST (b) model with noise coefficient $\mathbf{b} = [0, 1, 0]$. (c) Evolution of the mean modal energies, mean total energy, and mean helicity for the EST model with the strong noise coefficient $\mathbf{b} = [0, 10, 0]$</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>Fig. 5 Evolution of mean modal energies for a very large number of realisations for the EST (a) and HST (b) models. The mean is computed from 20,000 particles in the full noise case</figcaption></figure>

### 3.1.5 Data Assimilation

Using the two stochastic models in the full noise case described above, we perform the data assimilation tests using the following framework:

The signal process (the truth) is given by the deterministic triad model. The observations are the modal energies of the deterministic model, observed every 10 time units, and perturbed by noise of the form

$$ \eta \sim \mathcal{N}(\mathbf{0}, \mathbf{C}), \qquad (55) $$

where the covariance matrix $\mathbf{C} \in \mathbb{R}^3$ is chosen to be the diagonal matrix $\mathbf{C} = \text{diag}(0.005^2, 0.05^2, 0.05^2)$.

We use the sequential importance resampling (SIR) particle filter to assimilate the periodically observed signal process under the influence of observation noise. The particles evolve according to the stochastic triad models. Figures 6 and 7 show the results of filtering the ensemble of $n = 100$ particles of the EST and HST triad models, respectively. The ensembles are assessed in terms of the bias and RMSE statistics. We analyse the comparison details below:
<figure><img src="image_51.png" /><figcaption>(a) EST filter k</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>(b) EST filter p</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>(c) EST filter q</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(d) EST bias k</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>(e) EST bias p</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>(f) EST bias q</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>(g) EST RMSE k</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>(h) EST RMSE p</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>(i) EST RMSE q</figcaption></figure>

**Fig. 6** Filtering experiment for EST model using SIR particle filter. (a–c) Ensemble evolution for 100 particles in mode k (a, red), p (b, blue), and q (c, green). The signal (grey) is the deterministic model and the observations (black stars) are noisy and taken and assimilated every 10 time units. (d–f) Bias of the filtering ensemble. (g–h) RMSE of the filtering ensemble

## Mode k

This is the least energetic of all the modes (hence the reason why we observe it with the least amount of measurement noise). The cloud of particles is well placed around the truth even with the small sample. The bias remains small for both the HST and the EST versions and it reduces significantly when observations are assimilated more frequently (see Fig. 9 in Appendix 3) as well as when we use a large number (500) of particles (see Fig. 10 in Appendix 3). The RMSE remains small in all cases and decreases (though not substantially) when the DA step is small.

## Modes p and q

These are the two energetic modes of the system. We used here a measurement noise that is one order of magnitude larger. Despite this, the results remain equally good. The cloud of particles provide a good envelope for the truth at all times. This validates the choice of the stochasticity: the uncertainty is properly modelled.
<figure><img src="image_60.png" /><figcaption>**Fig. 7** Filtering experiment for HST model using SIR particle filter. (a–c) Ensemble evolution for 100 particles in mode k (a, red), p (b, blue), and q (c, green). The signal (grey) is the deterministic model and the observations (black stars) are noisy and taken and assimilated every 10 time units. (d–f) Bias of the filtering ensemble. (g–h) RMSE of the filtering ensemble</figcaption></figure>

For both models the bias can become very large, reaching 30% of the size of the oscillations for the HST model and 25% of the size of the oscillations for the EST model. As expected, it is drastically reduced when observations are assimilated more frequently. The RMSE for mode p is also large but substantially smaller for mode q. The addition of more intermediate DA steps or more particles has a less pronounced effect for the q mode.

**Remark 2** We record the Effective Sample Size (ESS) for a typical run (for both EST and HST) in Fig. 11. As usual, the ESS is computed just before the application of the resampling procedure. The ESS is seen to decay dramatically from 100 down to single digits numbers in most instances in time.

**Remark 3** We record the results over 10 independent runs of the filtering experiment for the EST model with 500 ensemble members in Fig. 12. More precisely, in graphs 12a, b and c each, we plot the mean across the 10 independent filtering runs together with the evolution of the signal and the individual ensemble means for each mode. The mean bias as well as the envelope obtained from the independent runs are shown in graphs 12d, e and f. The same is shown for the RMSE in graphs 12g,
h and i. Compared with a single run of the same experiment reported in Fig. 10 we observe the approximations are now near perfect (the statistical error has been drastically reduced).

# 4 Conclusions

The introduction of stochasticity into the deterministic triad models leads to two new stochastic models. Stochasticity is introduced in a principled way (rather than ad-hoc). It starts with a full scale fluid dynamic model which is randomly perturbed. At the full scale, the stochastic parametrisation models the small-scale effects in fluid dynamics modelling. In particular, it efficiently captures the high-frequency small-scale dynamics and correctly correlates it with the slow, large-scale fluid motion. In addition, it is constrained to conserve either the helicity or the kinetic energy of the system. This inspires two different stochastic triad models of Euler type which we compare using data assimilation procedures based on particle filtering. The methodology we employ can be used as a benchmark when analysing new types of stochastic parametrisations: ours is the first study that assesses the efficiency of stochastic parametrisations from a data assimilation perspective.

The introduction of stochasticity ensures that the correct spread (one that preserves the physical properties of the system) is introduced in the ensemble of particles. In its absence, the particle filter degenerates quite rapidly: after a few DA steps, a single particle survives the culling procedure which does offer a good approximation of the truth. A purely deterministic transition kernel does not work, generating a rapid degeneracy of the particle filter.

The two stochastic systems (one which preserves helicity and the other one which preserves energy) are analysed using a standard particle filter. There is no need for additional procedures (such as tempering, nudging, or jittering). They perform equally well from the viewpoint of DA: both the RMSE and the bias are drastically reduced and stabilised when the noise is carefully calibrated. The two different stochastic kernels require different noise calibrations in order to perform well in similar data assimilation scenarios. This is somehow expected, given that the underlying stochastic parametrisations preserve different physical quantities.

# Appendix 1: Notation and Basic Identities

## Notation

In this work we use the following notation. We write $\mathbb{Z}$, $\mathbb{R}$ and $\mathbb{C}$ for the sets of integers, real numbers and complex numbers, respectively. Boldface letters denote three-dimensional complex vectors. For two complex vectors $\mathbf{a}$ and $\mathbf{b}$ with
components $a_j$ and $b_j$, their dot product is denoted by

$$ \mathbf{a} \cdot \mathbf{b} = \sum_{j=1}^{3} a_j b_j = a_j b_j \in \mathbb{C}. $$

This paper uses the Einstein convention of summing over repeated indices. The norm of the complex vector **a** is defined as

$$ |\mathbf{a}| = \sqrt{\mathbf{a} \cdot \mathbf{a}^*} = \sqrt{a_j a_j^*} \ge 0, $$

with the superscript symbol * denoting complex conjugation. Further, the cross product of two vectors **a** and **b** is given by

$$ \mathbf{a} \times \mathbf{b} = (a_2 b_3 - a_3 b_2, a_3 b_1 - a_1 b_3, a_1 b_2 - a_2 b_1) \in \mathbb{C}^3. $$

The gradient of a scalar field $\phi : \mathcal{D} \subseteq \mathbb{R}^3 \to \mathbb{C}$ at a point $\mathbf{x} \in \mathcal{D}$ is denoted by

$$ \nabla \phi(\mathbf{x}) = (\partial_1 \phi(\mathbf{x}), \partial_2 \phi(\mathbf{x}), \partial_3 \phi(\mathbf{x})) \in \mathbb{C}^3. $$

The divergence of a vector field $\psi : \mathcal{D} \subseteq \mathbb{R}^3 \to \mathbb{C}^3$ with components $\psi_j$ at a point $\mathbf{x} \in \mathcal{D}$ is defined as

$$ \nabla \cdot \psi(\mathbf{x}) = \partial_1 \psi_1(\mathbf{x}) + \partial_2 \psi_2(\mathbf{x}) + \partial_3 \psi_3(\mathbf{x}) = \partial_j \psi_j(\mathbf{x}) \in \mathbb{C} $$

and the curl of $\psi$ at $\mathbf{x}$ is given by

$$ \nabla \times \psi(\mathbf{x}) = (\partial_2 \psi_3(\mathbf{x}) - \partial_3 \psi_2(\mathbf{x}), \partial_3 \psi_1(\mathbf{x}) - \partial_1 \psi_3(\mathbf{x}), \partial_1 \psi_2(\mathbf{x}) - \partial_2 \psi_1(\mathbf{x})) \in \mathbb{C}^3. $$

## Vector Identities

For three vectors **a**, **b** and **c**, we have the following algebraic vector identities

$$ \mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) = \mathbf{b} \cdot (\mathbf{c} \times \mathbf{a}) = \mathbf{c} \cdot (\mathbf{a} \times \mathbf{b}), $$

$$ \mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = (\mathbf{a} \cdot \mathbf{c})\mathbf{b} - (\mathbf{a} \cdot \mathbf{b})\mathbf{c}, $$

$$ (\mathbf{a} \times \mathbf{b}) \cdot (\mathbf{c} \times \mathbf{d}) = (\mathbf{a} \cdot \mathbf{c})(\mathbf{b} \cdot \mathbf{d}) - (\mathbf{b} \cdot \mathbf{c})(\mathbf{a} \cdot \mathbf{d}), $$

$$ \mathbf{a} \times \mathbf{a} = \mathbf{0}. $$

Moreover, we have the vector calculus identity

$$ (\mathbf{a} \cdot \nabla)\mathbf{b} + b_j \nabla a_j = -\mathbf{a} \times (\nabla \times \mathbf{b}) + \nabla(\mathbf{a} \cdot \mathbf{b}). \quad (56) $$
# Appendix 2: Derivation of Triad Models

## Deterministic Euler

We compute the projection of the terms corresponding to the deterministic Euler vorticity equation in curl form:

$$ \partial_t \mathbf{v} - \mathbf{v} \times (\nabla \times \mathbf{v}) + \frac{1}{2} \nabla |\mathbf{v}|^2 = -\nabla p. \quad (57) $$

onto the helical basis. For the time-derivative we get

$$ \begin{aligned} \int_{\mathcal{D}} \partial_t \mathbf{v}(\mathbf{x}, t) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{x}} \, \mathrm{d}\mathbf{x} &= \sum_{\mathbf{p}} \sum_{s_p} \partial_t a_{s_p}(\mathbf{p}, t) \mathbf{h}_{s_p}(\mathbf{p}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \int_{\mathcal{D}} e^{i(\mathbf{p}-\mathbf{k})\cdot\mathbf{x}} \, \mathrm{d}\mathbf{x} \\ &= L^3 \sum_{s_p} \partial_t a_{s_p}(\mathbf{k}, t) \mathbf{h}_{s_p}(\mathbf{k}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \\ &= L^3 \partial_t a_{s_k}(\mathbf{k}, t) \mathbf{h}_{s_k}(\mathbf{k}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \\ &= 2L^3 \partial_t a_{s_k}(\mathbf{k}, t). \end{aligned} $$

The vorticity term gives

$$ \begin{aligned} & \int_{\mathcal{D}} (\mathbf{v}(\mathbf{x}, t) \times \boldsymbol{\omega}(\mathbf{x}, t)) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{x}} \, \mathrm{d}\mathbf{x} \\ &= \sum_{\mathbf{p} \cdot \mathbf{q}} \sum_{s_p, s_q} a_{s_p}(\mathbf{p}, t) s_q |\mathbf{q}| a_{s_q}(\mathbf{q}, t) \mathbf{h}_{s_p}(\mathbf{p}) \times \mathbf{h}_{s_q}(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \int_{\mathcal{D}} e^{i(\mathbf{p}+\mathbf{q}-\mathbf{k})\cdot\mathbf{x}} \, \mathrm{d}\mathbf{x} \\ &= L^3 \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} a_{s_p}^*(\mathbf{p}, t) s_q |\mathbf{q}| a_{s_q}^*(\mathbf{q}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}). \end{aligned} \quad (58) $$

Note that we can write (58) in a form which is symmetric in $\mathbf{p}$ and $\mathbf{q}$ since, renaming the indices,

$$ \begin{aligned} & \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} a_{s_p}^*(\mathbf{p}, t) s_q |\mathbf{q}| a_{s_q}^*(\mathbf{q}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \\ &= \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} a_{s_q}^*(\mathbf{q}, t) s_p |\mathbf{p}| a_{s_p}^*(\mathbf{p}, t) \mathbf{h}_{s_q}^*(\mathbf{q}) \times \mathbf{h}_{s_p}^*(\mathbf{p}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \\ &= - \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} a_{s_q}^*(\mathbf{q}, t) s_p |\mathbf{p}| a_{s_p}^*(\mathbf{p}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}). \end{aligned} $$
Therefore,

$$ \int_D (\mathbf{v}(\mathbf{x}, t) \times \boldsymbol{\omega}(\mathbf{x}, t)) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{x}} \, d\mathbf{x} \quad (59) $$

$$ = \frac{L^3}{2} \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} (s_q |\mathbf{q}| - s_p |\mathbf{p}|) a_{s_p}^*(\mathbf{p}, t) a_{s_q}^*(\mathbf{q}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}). \quad (60) $$

Moreover, the gradient terms in (57) vanish upon expansion into helical modes. Thus, the Euler equation (57) in helical basis becomes

$$ \partial_t a_{s_k}(\mathbf{k}, t) = -\frac{1}{4} \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} (s_p |\mathbf{p}| - s_q |\mathbf{q}|) a_{s_p}^*(\mathbf{p}, t) a_{s_q}^*(\mathbf{q}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}). $$

## SALT Euler

We expand the 3D SALT Navier-Stokes equation (23) using (6) and (37). Assume $b_s(-\mathbf{p}) = b_s^*(\mathbf{p})$.

$$ \int_D \mathrm{d}\mathbf{v}(\mathbf{x}, t) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{x}} \, d\mathbf{x} = \sum_{\mathbf{p}} \sum_{s_p} \mathrm{d}a_{s_p}(\mathbf{p}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \int_D e^{i(\mathbf{p}-\mathbf{k}) \cdot \mathbf{x}} \, d\mathbf{x} \\ = 2L^3 \mathrm{d}a_{s_k}(\mathbf{k}, t) $$

And

$$ \begin{aligned} & \int_D \mathrm{d}\mathbf{x}_t(\mathbf{x}, t) \times \boldsymbol{\omega}(\mathbf{x}, t) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{x}} \, d\mathbf{x} \\ &= \sum_{\mathbf{p}, \mathbf{q}} \sum_{s_p, s_q} [a_{s_p}(\mathbf{p}, t) \mathrm{d}t + b_{s_p}(\mathbf{p}) \circ \mathrm{d}W_t] s_q |\mathbf{q}| a_{s_q}(\mathbf{q}, t) \mathbf{h}_{s_p}(\mathbf{p}) \\ & \quad \times \mathbf{h}_{s_q}(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \int_D e^{i(\mathbf{p}+\mathbf{q}-\mathbf{k}) \cdot \mathbf{x}} \, d\mathbf{x} \\ &= L^3 \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} [a_{s_p}^*(\mathbf{p}, t) \mathrm{d}t + b_{s_p}^*(\mathbf{p}) \circ \mathrm{d}W_t] s_q |\mathbf{q}| a_{s_q}^*(\mathbf{q}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \\ & \quad \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}). \end{aligned} $$

Renaming the indices, we can write
$$
\begin{align*}
& L^3 \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} [a_{s_p}^*(\mathbf{p}, t)dt + b_{s_p}^*(\mathbf{p}) \circ dW_t]_{s_q} |\mathbf{q}| a_{s_q}^*(\mathbf{q}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \\
&= L^3 \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} [a_{s_q}^*(\mathbf{q}, t)dt + b_{s_q}^*(\mathbf{q}) \circ dW_t]_{s_p} |\mathbf{p}| a_{s_p}^*(\mathbf{p}, t) \mathbf{h}_{s_q}^*(\mathbf{q}) \times \mathbf{h}_{s_p}^*(\mathbf{p}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) \\
&= -L^3 \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} [a_{s_q}^*(\mathbf{q}, t)dt + b_{s_q}^*(\mathbf{q}) \circ dW_t]_{s_p} |\mathbf{p}| a_{s_p}^*(\mathbf{p}, t) \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}).
\end{align*}
$$

Thus, we arrive at

$$
\int_D d\mathbf{x}_t(\mathbf{x}, t) \times \omega(\mathbf{x}, t) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{x}} d\mathbf{x} \\
= \frac{L^3}{2} \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} \left[ (s_q |\mathbf{q}| b_{s_p}^*(\mathbf{p}) \circ dW_t) a_{s_q}^*(\mathbf{q}, t) - s_p |\mathbf{p}| b_{s_q}^*(\mathbf{q}) \circ dW_t) a_{s_p}^*(\mathbf{p}, t) \right. \\
\left. + (a_{s_p}^*(\mathbf{p}, t) dt s_q |\mathbf{q}| a_{s_q}^*(\mathbf{q}, t) - a_{s_q}^*(\mathbf{q}, t) dt s_p |\mathbf{p}| a_{s_p}^*(\mathbf{p}, t)) \right] \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}).
$$

Therefore,

$$
\begin{align*}
\mathrm{d}a_{s_k}(\mathbf{k}, t) ={}& \frac{1}{4} \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} \left[ (s_q |\mathbf{q}| b_{s_p}^*(\mathbf{p}) \circ dW_t) a_{s_q}^*(\mathbf{q}, t) - s_p |\mathbf{p}| b_{s_q}^*(\mathbf{q}) \circ dW_t) a_{s_p}^*(\mathbf{p}, t) \right. \\
& \left. \qquad + (a_{s_p}^*(\mathbf{p}, t) \mathrm{d}t s_q |\mathbf{q}| a_{s_q}^*(\mathbf{q}, t) - a_{s_q}^*(\mathbf{q}, t) \mathrm{d}t s_p |\mathbf{p}| a_{s_p}^*(\mathbf{p}, t)) \right] \mathbf{h}_{s_p}^*(\mathbf{p}) \times \mathbf{h}_{s_q}^*(\mathbf{q}) \cdot \mathbf{h}_{s_k}^*(\mathbf{k}).
\end{align*}
$$

LU Euler

Written in terms of the SALT model, the LU model is

$$
\mathrm{d}\mathbf{v} + \mathrm{d}\mathbf{x}_t \cdot \nabla \mathbf{v} + \mathbf{v}^j \nabla \mathbf{x}_t^j - \mathbf{v}^j \nabla (\xi \circ dW_t)^j = -\nabla \mathrm{d}p.
$$

Expanding the additional term gives

$$
\nabla (\xi \circ dW_t)^j = \nabla \sum_{\mathbf{q}} (b_{\pm}(\mathbf{q})\mathbf{h}_{\pm}(\mathbf{q})e^{i\mathbf{q} \cdot \mathbf{x}} \circ dW_t)^j = \sum_{\mathbf{q}} i\mathbf{q}b_{\pm}(\mathbf{q})\mathbf{h}_{\pm}^j(\mathbf{q})e^{i\mathbf{q} \cdot \mathbf{x}} \circ dW_t.
$$
Thus we get

$$ \mathbf{v}^j \nabla (\xi \circ dW_t)^j = \sum_{\mathbf{p}, \mathbf{q}} \sum_{s_p, s_q} a_{s_p}(\mathbf{p}, t) \mathbf{h}_{s_p}^j(\mathbf{p}) i \mathbf{q} b_{s_q}(\mathbf{q}) \mathbf{h}_{s_q}^j(\mathbf{q}) e^{i(\mathbf{p}+\mathbf{q})\cdot\mathbf{x}} \circ dW_t. $$

Now projecting and renaming, we have

$$ \begin{aligned} \int_D \mathbf{v}^j \nabla (\xi \circ dW_t)^j \cdot a_{s_k}^*(\mathbf{k}, t) \mathbf{h}_{s_k}^*(\mathbf{k}) e^{-i\mathbf{k} \cdot \mathbf{x}} \, d\mathbf{x} &= \sum_{\mathbf{p}, \mathbf{q}} \sum_{s_p, s_q} a_{s_k}^*(\mathbf{k}, t) a_{s_p}(\mathbf{p}, t) [b_{s_q}(\mathbf{q}) \circ dW_t] (i\mathbf{q} \cdot \mathbf{h}_{s_k}^*(\mathbf{k})) (\mathbf{h}_{s_p}(\mathbf{p}) \cdot \mathbf{h}_{s_q}(\mathbf{q})) \int_D e^{(\mathbf{p}+\mathbf{q}-\mathbf{k})\cdot\mathbf{x}} \, d\mathbf{x} \\ &= L^3 \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} a_{s_k}^*(\mathbf{k}, t) a_{s_p}^*(\mathbf{p}, t) [b_{s_q}^*(\mathbf{q}) \circ dW_t] (-i\mathbf{q} \cdot \mathbf{h}_{s_k}^*(\mathbf{k})) (\mathbf{h}_{s_p}^*(\mathbf{p}) \cdot \mathbf{h}_{s_q}^*(\mathbf{q})) \\ &= L^3 \sum_{\mathbf{p}+\mathbf{q}+\mathbf{k}=0} \sum_{s_p, s_q} (f_k^{pq})^* a_{s_k}^*(\mathbf{k}, t) (a_{s_p}^*(\mathbf{p}, t) [b_{s_q}^*(\mathbf{q}) \circ dW_t] + a_{s_q}^*(\mathbf{q}, t) [b_{s_p}^*(\mathbf{p}) \circ dW_t]) \end{aligned} $$

with

$$ f_k^{pq} = (-i(\mathbf{p} + \mathbf{q}) \cdot \mathbf{h}_{s_k}(\mathbf{k})) (\mathbf{h}_{s_p}(\mathbf{p}) \cdot \mathbf{h}_{s_q}(\mathbf{q})). $$

Note that, due to the triad condition $\mathbf{p} + \mathbf{q} + \mathbf{k} = 0$,

$$ f_k^{pq} = f_p^{kq} = f_q^{kp} = 0 $$

so that the difference term between SALT and LU vanishes in the helical projection and the two projected models coincide.

## Appendix 3: Supplementary Numerics

### Calibration of the Noise Amplitude

To calibrate the noise amplitude for the data assimilation experiments, we rely on two forecast verification metrics. The rank histogram (or Talagrand histogram) and the continuous ranked probability score (CRPS). We evaluated these for both models on 64 different noise amplitude vectors. The metrics were recorded by running the data assimilation/particle filtering experiment described in the main text for the noise amplitude vectors $\mathbf{b} = [b_k, b_p, b_q]$ resulting from all possible combinations of $b_k \in \{0.05, 0.1, 0.2, 0.5\}$, $b_p \in \{0.025, 0.05, 0.1, 0.2\}$ and $b_q \in \{0.01, 0.02, 0.04, 0.1\}$. We present the mean CRPS scores for the top 5 tested noise
<table><caption>**Table 1** CRPS scores (lower is better) for the five best—in terms of overall mean CRPS score—tested noise amplitude vectors **b** (first column). The mean CRPS score over the three modal energies is shown for the EST (second column) and HST (third column) models. The overall mean CRPS score for the respective noise amplitude is shown in the *Mean* column (fourth column). Finally, we provide references to the figures showing the associated rank histograms (last column)</caption><thead><tr><th>b</th><th>EST</th><th>HST</th><th>Mean</th><th>Histograms</th></tr></thead><tbody><tr><td>[0.05, 0.025, 0.01]</td><td>0.0282</td><td>0.0394</td><td>0.0338</td><td>Figure 8a + f</td></tr><tr><td>[0.05, 0.025, 0.02]</td><td>0.0338</td><td>0.0392</td><td>0.0365</td><td>Figure 8b + g</td></tr><tr><td>[0.10, 0.025, 0.02]</td><td>0.0352</td><td>0.0385</td><td>0.0368</td><td>Figure 8c + h</td></tr><tr><td>[0.10, 0.05, 0.02]</td><td>0.0393</td><td>0.0356</td><td>0.0375</td><td>Figure 8d + i</td></tr><tr><td>[0.10, 0.05, 0.01]</td><td>0.0363</td><td>0.0409</td><td>0.0386</td><td>Figure 8e + j</td></tr></tbody></table>

<figure><img src="image_61.png" /><figcaption>**Fig. 8** Rank histograms for the 5 best noise amplitude vectors in terms of CRPS score (see Table 1) for the EST and HST models. Each individual graph shows the rank histogram of an ensemble of 15 particles run up to a final time of 1400, with data assimilation performed every 10 time units. The top histogram in each subfigure represents the ensemble for mode **k**, the middle histogram represents the ensemble for mode **p** and the bottom one for mode **q**</figcaption></figure>

amplitude vectors in Table 1. The mean is taken across both models. Based on this calibration, we chose the case **b** = [0.10, 0.05, 0.01] for the data assimilation experiment, as the other ones have inferior rank histograms, so we achieve a good balance between the visual judgment of rank histograms and the CRPS score (Fig. 8).

## Data Assimilation Verification

See Figs. 9, 10, 11, and 12.
<figure><img src="image_62.png" /><figcaption>(a) EST filter k</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>(b) EST filter p</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>(c) EST filter q</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>(d) EST bias k</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>(e) EST bias p</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>(f) EST bias q</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>(g) EST RMSE k</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>(h) EST RMSE p</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>(i) EST RMSE q</figcaption></figure>

**Fig. 9** Filtering experiment for EST model using SIR particle filter with a small data assimilation interval of 5 time units. (a–c) Ensemble evolution for 100 particles in mode k (a, red), p (b, blue), and q (c, green). The signal (grey) is the deterministic model and the observations (black stars) are noisy and taken and assimilated every 5 time units. (d–f) Bias of the filtering ensemble. (g–h) RMSE of the filtering ensemble
<figure><img src="image_71.png" /><figcaption>(a) EST filter k</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>(b) EST filter p</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>(c) EST filter q</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>(d) EST bias k</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>(e) EST bias p</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>(f) EST bias q</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>(g) EST RMSE k</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>(h) EST RMSE p</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>(i) EST RMSE q</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>(a) ESS for EST experiment</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>(b) ESS for HST experiment</figcaption></figure>

<caption>**Fig. 11** Typical ESS for the filtering experiments. See Remark 2 in the main text. (a) ESS for EST experiment. (b) ESS for HST experiment</caption>
<figure><img src="image_82.png" /><figcaption>**Fig. 12** Statistics of the mean ensemble over 10 independent runs of the particle filter with 500 ensemble members for the EST model. See Remark 3 in the main text</figcaption></figure>

References

[1] G. Tadmor B. Noack, M. Morzinsky. *Reduced-Order Modelling for Flow Control*, volume 528. CISM Courses and Lectures Springer, 2010.

[2] Alan Bain and Dan Crisan. *Fundamentals of stochastic filtering*, volume 60 of *Stochastic Modelling and Applied Probability*. Springer, New York, 2009. Available from: https://doi.org/10.1007/978-0-387-76896-0.

[3] Carl M Bender, Darryl D Holm, and Daniel W Hook. Complexified dynamical systems. *Journal of Physics A: Mathematical and Theoretical*, 40(32):F793, 2007.

[4] Judith Berner, Ulrich Achatz, Lauriane Batte, Lisa Bengtsson, Alvaro De La Camara, Hannah M Christensen, Matteo Colangeli, Danielle RB Coleman, Daan Crommelin, Stamen I Dolaptchiev, et al. Stochastic parameterization: Toward a new view of weather and climate models. *Bulletin of the American Meteorological Society*, 98(3):565–588, 2017.

[5] L. Biferale. Shell models of energy cascade in turbulence. *Annual Review of Fluid Mechanics*, 35:441–468, 2003.

[6] L. Biferale. Shell models of energy cascade in turbulence. *Annual Rev. Fluid Mech.*, 35:441–468, 2003.

[7] Tomas Bohr, Mogens H Jensen, Giovanni Paladin, and Angelo Vulpiani. *Dynamical systems approach to turbulence*. CUP, 1998.
[8] B. Chapron, P. Dérion, E. Mémin, and V. Ressaguier. Large-scale flows under location uncertainty: a consistent stochastic framework. *Quarterly Journal of the Royal Meteorological Society*, 144(710):251–260, 2018. Available from: https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/qj.3198, arXiv:https://rmets.onlinelibrary.wiley.com/doi/pdf/10.1002/qj.3198, https://doi.org/10.1002/qj.3198.

[9] Qiaoning Chen, Shiyi Chen, and Gregory L Eyink. The joint cascade of energy and helicity in three-dimensional turbulence. *Physics of Fluids*, 15(2):361–374, 2003. Available from: https://doi.org/10.1063/1.1533070.

[10] Qiaoning Chen, Shiyi Chen, Gregory L Eyink, and Darryl D Holm. Intermittency in the joint cascade of energy and helicity. *Physical Review Letters*, 90(21):214503, 2003. Available from: https://doi.org/10.1103/PhysRevLett.90.214503.

[11] Colin J. Cotter, Dan Crisan, Darryl D. Holm, Wei Pan, and Igor Shevchenko. Numerically modelling stochastic lie transport in fluid dynamics, 2018. Available from: https://arxiv.org/abs/1801.09729, https://doi.org/10.48550/ARXIV.1801.09729.

[12] Alex DD Craik. *Wave interactions and fluid flows*. Cambridge Nonlinear Sciences, 1998.

[13] B. Geurts, D. Holm, and E. Luesink. Lyapunov exponents of two stochastic lorenz 63 systems. *Journal of statistical physics*, 179:1343–1365, 2017.

[14] Bernard J. Geurts, Darryl D. Holm, and Erwin Luesink. Lyapunov exponents of two stochastic Lorenz 63 systems. *Journal of Statistical Physics*, 179(5):1343–1365, Jun 2020. Available from: https://doi.org/10.1007/s10955-019-02457-3.

[15] E. B. Gledzer. Rotation and helicity effects in cascade models of turbulence. *Doklady Physics*, 53:216–220, 2008.

[16] Darryl D. Holm. Variational principles for stochastic fluid dynamics. *Proc. A.*, 471(2176):20140963, 19, 2015.

[17] Darryl D Holm and Ruiao Hu. Stochastic effects of waves on currents in the ocean mixed layer. *Journal of Mathematical Physics*, 62(7):073102, 2021.

[18] Darryl D Holm, Jerrold E Marsden, and Tudor S Ratiu. The Euler-Poincaré equations and semidirect products with applications to continuum theories. *Advances in Mathematics*, 137(1):1–81, 1998.

[19] Ioannis Karatzas and Steven Shreve. *Brownian motion and stochastic calculus*, volume 113. Springer Science & Business Media, 2012.

[20] Robert H Kraichnan. Small-scale structure of a scalar field convected by turbulence. *The Physics of Fluids*, 11(5):945–953, 1968.

[21] Long Li. *Stochastic modeling and numerical simulation of ocean dynamics*. PhD thesis, Université de Rennes 1, 2021.

[22] Alexander Lobbe. alobbe/stochastic-triads: Initial release, April 2023. Available from: https://doi.org/10.5281/zenodo.7845270.

[23] Etienne Mémin. Fluid flow dynamics under location uncertainty. *Geophys. Astrophys. Fluid Dyn.*, 108(2):119–146, 2014.

[24] E. A. Novikov. *Infinitely divisible distributions in turbulence*. Physical Review E, 1994.

[25] Valentin Ressaguier, Agustin M. Picard, Etienne Memin, and Bertrand Chapron. Quantifying truncation-related uncertainties in unsteady fluid dynamics reduced order models. *SIAM/ASA Journal on Uncertainty Quantification*, 9(3):1152–1183, 2021. Available from: https://doi.org/10.1137/19M1354819.

[26] Peter Jan van Leeuwen, Dan Crisan, Oana Lang, and Roland Potthast. Bayesian inference for fluid dynamics: A case study for the stochastic rotating shallow water model, 2021. arXiv:2112.15216.

[27] Dario Vincenzi and John D Gibbon. How close are shell models to the 3d Navier-Stokes equations? *Nonlinearity*, 34(8):5821, 2021.

[28] Fabian Waleffe. The nature of triad interactions in homogeneous turbulence. *Physics of Fluids A: Fluid Dynamics*, 4(2):350–363, 1992.

[29] Fabian Waleffe. Inertial transfers in the helical decomposition. *Physics of Fluids A: Fluid Dynamics*, 5(3):677–685, 1993. Available from: https://doi.org/10.1063/1.858651.
**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_83.png" /></figure>
# An Explicit Method to Determine Casimirs in 2D Geophysical Flows

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

Check for
updates

Erwin Luesink and Bernard Geurts

**Abstract** Conserved quantities in geophysical flows play an important role in the characterisation of geophysical dynamics and aid the development of structure-preserving numerical methods. A significant family of conserved quantities is formed by the Casimirs i.e., integral conservation laws that are in the kernel of the underlying Poisson bracket. The Casimirs hence determine the geometric structure of the geophysical fluid equations among which the enstrophy is well known. Often Casimirs are proposed on heuristic grounds and later verified to be part of the kernel of the Poisson bracket. In this work, we will explicitly construct Casimirs by rewriting the Poisson bracket in vorticity-divergence coordinates thereby providing explicit construction of Casimirs for 2D geophysical flow dynamics.

## 1 Introduction

Models of geophysical flows involve a fluid with a free interface under the influence of gravity on a rotating domain, often forced by variations in temperature, salinity and density. In this paper, we will consider the thermal rotating shallow water equations, which is a two dimensional model that includes all of the features above.

E. Luesink (✉)

Multiscale Modelling and Simulation, Department of Applied Mathematics, Faculty EEMCS,
University of Twente, Enschede, The Netherlands
e-mail: e.luesink@utwente.nl

B. Geurts

Multiscale Modelling and Simulation, Department of Applied Mathematics, Faculty EEMCS,
University of Twente, Enschede, The Netherlands

Multiscale Physics, Center for Computational Energy Research, Department of Applied Physics,
Eindhoven University of Technology, Eindhoven, The Netherlands

<footnote>This work was performed in the project SPRESTO (structure-preserving regularization and stochastic forcing for nonlinear hyperbolic PDEs), supported by a NWO TOP 1 grant.</footnote>
The variations of temperature and salinity are collected in a single scalar field that we call the buoyancy. The thermal rotating shallow water equations can be derived in several different ways from the rotating Euler equations with stratification. The best known method is to assume that the domain is shallow, allowing the replacement of the prognostic equation for the vertical velocity by a hydrostatic pressure relation from which the vertical velocity can be inferred. The result is the system of inviscid primitive equations. Upon assuming columnar motion, e.g., motivated by the Taylor-Proudman theorem, meaning that the horizontal velocity field does not depend on the vertical coordinate, one can vertically average the primitive equations to obtain the thermal rotating shallow water equations.

The focus in this contribution is with the preservation of key structures that characterise the governing equations and their dynamics. In geophysical fluid dynamics (GFD) models in 2D so-called conservation of Casimirs forms a major framework for the modelling and construction of methods for their numerical treatment. We present an explicit constructive method with which the Casimirs of a general class of GFD models can be explicitly computed. This provides an alternative to earlier work by Cotter and Holm [2013] that yields explicit Casimirs for 2D systems.

This entire sequence of assumptions and approximations can be performed at the level of the Lagrangian formulation of the models, as shown in Holm and Luesink [2021]. Doing so is helpful because the geometric structure associated with such a variational formulation implies various important dynamical quantities and conservation laws. For two-dimensional geophysical fluids, we find that the key dynamical quantities are the potential vorticity and the potential buoyancy, where the latter is the ratio of the buoyancy to the depth. Both quantities play a key role in cyclogenesis, as shown in Holm et al. [2021]. In the present work, we will show that the potential vorticity and potential buoyancy can also be used to give an alternative formulation of the fluid equations. This formulation is particularly useful for the identification of Casimirs, which are conserved integral quantities. The Casimirs in turn help interpret mechanisms in geophysical turbulence. For instance, in the incompressible setting without buoyancy effects in two dimensions, the entrophy is one of the Casimirs and plays a central role in the double cascade predicted by Kraichnan [1967]. From a geometric point of view, the Casimirs are the functionals that form the kernel of the Poisson bracket. Usually, the form of the Casimirs is assumed or guessed and then verified by checking whether the Poisson bracket indeed vanishes. We will use two subsequent changes of variables to provide a constructive derivation for the form of the Casimirs.

The organization of this paper is as follows. In Sect. 2 we present the main equations relevant to 2D geophysical flow and sketch the central importance of the retaining symmetries and conserved properties of the models. Among these properties, Casimirs are a major structure for the GFD models in 2D. Explicit methods to determine these Casimirs will be sketched in Sect. 3. Concluding remarks are collected in Sect. 4.
## 2 Geophysical Flows

In Holm and Luesink [2021] it is shown that starting from the Lagrangian for the rotating, stratified Euler equations in three dimensions, one can derive the Lagrangian for the thermal rotating shallow water equations by assuming hydrostatic pressure and applying vertical averaging. The benefit of performing this sequence of approximations at the level of the Lagrangian is that the geometric structure is not affected. This means that the thermal rotating shallow water equations can be formulated using a Lie-Poisson bracket. As was shown in Marsden and Weinstein [1983], there exist multiple equivalent formulations in terms of Lie-Poisson brackets of the two-dimensional Euler equations. Similar formulations can be derived for the compressible two-dimensional fluid models. This was first done in Holm and Long [1989]. The retained Lagrangian structure in 2D is basic to a possible systematic approach to determining the conserved quantities associated with the thermal rotating shallow water equations. This systematic approach will be elaborated on in this paper.

It is possible to derive the equations for geophysical fluid dynamics on arbitrary smooth manifolds starting from the Euler equations, as shown in Holm et al. [1998]. This level of generality requires tools of differential topology, such as the Lie derivative, the pullback and the pushforward, see for instance Abraham and Marsden [1978], Marsden and Ratiu [2013], Holm et al. [2009]. However, if the domain is a two-dimensional compact subset $\Omega$ embedded in $\mathbb{R}^3$ or $\mathbb{R}^2$ with Cartesian coordinates and appropriate boundary conditions, we can formulate the equations of motion also using vector proxies of exterior calculus. The computations are in arbitrary orthonormal coordinate systems. To this end, we start from the dimensionless Lagrangian for the thermal rotating shallow water (TRSW) equations, which are given by

$$L_{\text{TRSW}}(\mathbf{u}, \eta, b) = \int_{\Omega} \left( \frac{1}{2} |\mathbf{u}|^2 + \frac{1}{\text{Ro}} \mathbf{u} \cdot \mathbf{R} - \frac{1}{2 \text{Fr}^2} (1 + \varepsilon b)(\eta - 2h) \right) \eta \, d\mu. \quad (1)$$

In this Lagrangian, $\mathbf{u}$ is the velocity field, $\eta = \alpha\zeta + h$ is the total depth, $\zeta$ is the free surface elevation, $h$ is the bottom topography, $\mathbf{R}$ is the vector potential for the Coriolis parameter and $b$ is the buoyancy variable. Furthermore, the dimensionless numbers are the aspect ratio $\sigma$, the Rossby number $\text{Ro}$, the Froude number $\text{Fr}$, the wave amplitude $\alpha$ and the stratification parameter $\varepsilon$. The stratification parameter represents the importance of the buoyancy variable, which itself is of order one. The vector potential for the Coriolis parameter satisfies $\nabla^{\perp} \cdot \mathbf{R} = f(x, y)$, where $f(x, y)$ is the usual Coriolis parameter. The $\perp$ operator corresponds to a Hodge star operator, which for Cartesian coordinates is defined by $(x, y)^{\perp} = (-y, x)$. The volume form is given by $d\mu$, which in Cartesian coordinates is expressed as $d\mu = dx dy$. The Lagrangian is a functional on $\mathcal{X}(\Omega) \times V^*(\Omega)$, the product space of the space of vector fields $\mathcal{X}(\Omega)$ and the space of advected quantities $V^*(\Omega)$. That is, $\mathbf{u} \in \mathcal{X}(\Omega)$
and $\eta, b \in V^*(\Omega)$. More information and details on these spaces can be found in Luesink [2021]. The thermal rotating shallow water equations associated with the Lagrangian (1) are given in the advective formulation by

$$ 
\begin{aligned}
& \frac{\partial}{\partial t} \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \frac{1}{\text{Ro}} f \mathbf{u}^{\perp} = -\frac{\alpha}{\text{Fr}^2} \nabla((1 + \text{sh})\zeta) + \frac{\zeta}{2 \text{Fr}^2} (\alpha\zeta - h) \nabla b, \\
& \frac{\partial}{\partial t} \eta + \nabla \cdot (\eta \mathbf{u}) = 0, \\
& \frac{\partial}{\partial t} b + (\mathbf{u} \cdot \nabla) b = 0.
\end{aligned}
\quad (2) $$

Equations (2) describe a compressible fluid with thermal effects in a rotating frame in two dimensions. Note that the equations take the same form on the sphere. The bold symbols denote vector valued quantities, which assumes the existence of a basis. Since we are working with two dimensional manifolds, such a basis is available. We will now recall an important vector calculus identity that is central to the manipulations that we will perform. This identity follows from differential topology, as shown in Holm et al. [2009], where vector fields and differential forms in arbitrary coordinate systems can be defined. The vector-valued coefficients of these objects can be identified when the underlying space is Euclidean with the standard inner product. This means that $\mathbf{u}$ is simultaneously the coefficient of a vector field and a differential 1-form. This can lead to confusion, so in what follows $\mathbf{u}$ denotes a vector field in arbitrary coordinates and $\mathbf{v}$ denotes a differential 1-form. The Lie derivative can be defined as the derivative of the pullback

$$ \mathcal{L}_{\mathbf{u}} v = \frac{d}{dt} \Big|_{t=0} \phi_t^*(v), \quad (3) $$

where $\phi$ is the flow associated with the vector field $\mathbf{u}$. The asterisk on $\phi$ means that $v$ is pulled back along $\phi$. Here, $\phi$ is the solution to the differential equation $d\phi/dt = \mathbf{u}(\phi, t)$ with arbitrary initial condition, i.e., representing a trajectory of a test-particle in the flow field $\mathbf{u}$. The arbitrariness is what allows Lagrangian formulations of fluid dynamics to be related to Eulerian formulations. The Lie derivative can also be defined by means of Cartan's formula

$$ \mathcal{L}_{\mathbf{u}} v = i_{\mathbf{u}} dv + d(i_{\mathbf{u}} v), \quad (4) $$

where $i_{\mathbf{u}} v$ denotes the interior product of $\mathbf{u}$ with $v$ and $d$ is the exterior derivative, see e.g. Abraham and Marsden [1978], Marsden and Ratiu [2013], Holm et al. [2009]. Setting the two definitions (3) and (4) of the Lie derivative equal to one another in $\mathbb{R}^2$ provides the following vector calculus identity

$$ (\mathbf{u} \cdot \nabla)\mathbf{v} + (\nabla\mathbf{u}) \cdot \mathbf{v} = (\nabla^{\perp} \cdot \mathbf{v})\mathbf{u}^{\perp} + \nabla(\mathbf{u} \cdot \mathbf{v}), \quad (5) $$
where the left hand side results from the pullback definition (3) and the right hand side follows from the Cartan formula (4).

Using (5), an alternative formulation to (2) of the thermal rotating shallow water equations is given by

$$
\begin{align}
\frac{\partial}{\partial t} \boldsymbol{u} + \left( \omega + \frac{1}{\mathrm{Ro}} f \right) \boldsymbol{u}^{\perp} &= -\frac{\alpha}{\mathrm{Fr}^2} \nabla \left( (1 + \mathrm{s}b)\zeta \right) - \frac{1}{2} \nabla |\boldsymbol{u}|^2 + \frac{\mathrm{s}}{2 \mathrm{Fr}^2} (\alpha \zeta - h) \nabla b, \nonumber \\
\frac{\partial}{\partial t} \eta + \nabla \cdot (\eta \boldsymbol{u}) &= 0, \nonumber \\
\frac{\partial}{\partial t} b + (\boldsymbol{u} \cdot \nabla) b &= 0, \tag{6}
\end{align}
$$

where $\omega = \nabla^{\perp} \cdot \mathbf{u}$ is the scalar vorticity. This formulation is often called the vector invariant formulation of fluid dynamics. There are several situations where the vector invariant formulation has advantages over the advective formulation. The vector invariant form offers an easier basis for numerical discretisations because it hides the nonlinearity. Additionally, on the right hand side of the velocity equation in (6), one recognises the gradient of the Bernoulli function $B = \frac{1}{2}|\mathbf{u}|^2 + \frac{\alpha}{\mathrm{Fr}^2}(1+\mathrm{s}b)\zeta$, further explaining the physical mechanisms responsible for the flow of fluid.

Both (2) and (6) can be formulated as a Hamiltonian system with a Lie-Poisson bracket. This is a crucial property of these models as it prepares for a clear formulation of conserved quantities, i.e., the desired Casimir s. By Legendre transforming the Lagrangian (1), we obtain the Hamiltonian associated with the thermal rotating shallow water equations. The Legendre transform is given by

$$
H(\mathbf{m}, \eta, b) = \langle \mathbf{m}, \mathbf{u} \rangle - L(\mathbf{u}, \eta, b), \quad (7)
$$

where **m** = *δL*/(*δu*) is the momentum variable in terms of the functional derivative
of the Lagrangian. The Legendre transform only relates the velocity and the
momentum and does not affect the advected quantities. The Hamiltonian is therefore
a functional on the space *X*\*(Ω) × *V*\*(Ω), where *X*\*(Ω) is the dual space of *X*(Ω)
with respect to *L*²-pairing on Ω. For the thermal rotating shallow water model, the
momentum is given by **m** = *η*(**u** + **R**) and the Hamiltonian is

$$
H_{\text{TRSW}}(\mathbf{m}, \eta, b) = \int_{\Omega} \left( \frac{1}{2\eta^2} |\mathbf{m}|^2 + \frac{1}{2\text{Fr}^2} (1 + \text{s}b)(\eta - 2h) \right) \eta \, d\mu. \quad (8)
$$
The Poisson bracket that will allow us to formulate the equations of motion has the following general form, see e.g. Holm et al. [2009],

$$ \{F, G\} = \int_{\Omega} \begin{pmatrix} \delta F / \delta \mathbf{m} \\ \delta F / \delta \eta \\ \delta F / \delta b \end{pmatrix} \cdot \mathcal{J}(\mathbf{m}, \eta, b) \begin{pmatrix} \delta G / \delta \mathbf{m} \\ \delta G / \delta \eta \\ \delta G / \delta b \end{pmatrix} d\mu, \quad (9) $$

where $F$ and $G$ are two arbitrary functionals on the space $\mathcal{X}^*(\Omega) \times V^*$. The matrix $\mathcal{J}$ is operator valued and depends on the state. $\mathcal{J}$ can be explicitly derived using reduction theorems, following Marsden and Weinstein [1974], Holm et al. [1998]. We speak of Hamiltonian dynamics if the evolution of a functional $F$ is expressed by

$$ \frac{d}{dt} F = -\{F, H\}, \quad (10) $$

where $H$ is the Hamiltonian of the system. We obtain the equations of motion (2) by choosing $G = H$ and using $\mathcal{J}(\mathbf{m}, \eta, b)$ as in Holm et al. [2021] to find

$$ \frac{\partial}{\partial t} \begin{pmatrix} m_i \\ \eta \\ b \end{pmatrix} = \underbrace{\begin{pmatrix} m_j \partial_i + \partial_j m_i & \eta \partial_i - b_{,i} \\ \partial_j \eta & 0 \\ b_{,j} & 0 \end{pmatrix}}_{=\mathcal{J}(\mathbf{m}, \eta, b)} \begin{pmatrix} \delta H_{\text{TRSW}} / \delta m_j \\ \delta H_{\text{TRSW}} / \delta \eta \\ \delta H_{\text{TRSW}} / \delta b \end{pmatrix} \quad (11) $$

and reproduces Eqs. (2) after using the definition of the momentum, rewriting the momentum equation using the continuity equation and applying the Lie derivative identity (5) to the vector potential of the Coriolis parameter. We have used Einstein's summation convention of summing over repeated indices and the notation $b_{,i}$ refers to $i$-th component of the gradient of $b$ (the comma denotes a spatial derivative). Note that the dependence on the state $(\mathbf{m}, \eta, b)$ is linear. This means that the Poisson bracket is in fact a Lie-Poisson bracket. The kernel of a Lie-Poisson bracket is nontrivial and is precisely the kernel of the linear operator $\mathcal{J}$. The kernel is spanned by functionals known as Casimir and the goal of the remainder of this work is to derive these Casimir functionals explicitly.

In Holm and Long [1989], the momentum bracket (11) is transformed by means of a change of variables to the rotational form. The change of variables is $(m_1, m_2, \eta, b) \mapsto (u, v, \eta, b)$, where $u = \frac{m_1}{\eta}$ and $v = \frac{m_2}{\eta}$, which is invertible as long as $0 < \eta < \infty$. This change of variables is not the same as just formulating the Hamiltonian and the Lie-Poisson bracket in terms of velocity because of the role of the rotating frame. Performing the transformation of Holm and Long [1989], we
obtain the following formulation

$$
\frac{\partial}{\partial t} \begin{pmatrix} u \\ v \\ \eta \\ b \end{pmatrix} = \underbrace{\begin{pmatrix} 0 & -q & \partial_x - \frac{b_x}{\eta} \\ q & 0 & \partial_y - \frac{b_y}{\eta} \\ \partial_x & \partial_y & 0 \\ \frac{b_x}{\eta} & \frac{b_y}{\eta} & 0 \end{pmatrix}}_{=\mathcal{J}(u,v,\eta,b)} \begin{pmatrix} \delta H_{TRSW} / \delta u = \eta u \\ \delta H_{TRSW} / \delta v = \eta v \\ \delta H_{TRSW} / \delta \eta = B \\ \delta H_{TRSW} / \delta b = T \end{pmatrix}, \quad (12)
$$

where $q = \frac{1}{\eta}(\omega + \frac{1}{\text{Ro}}f)$ is the potential vorticity, $B = \frac{1}{2}|\mathbf{u}|^2 + \frac{1}{2\text{Fr}^2}(1+s\mathbf{b})(\eta - 2h)$ is the Bernoulli function and $T = \frac{s}{2\text{Fr}^2}(\eta^2 - 2\eta h)$. This Hamiltonian formulation corresponds to (6) and extends the bracket of Holm and Long [1989] to include thermal variations. This bracket is linear and it is a simple exercise to show that it is skew-symmetric with respect to the $L^2$-pairing. To prove that it satisfies the Jacobi identity, we repeat the argument of Holm and Long [1989] and state that the bracket in (12) is an invertible transformation of variables in the Lie-Poisson bracket (11).

We now wish to express the bracket (12) in another set of variables to achieve a formulation in independent scalar variables. This is accomplished using the transformation $(u, v, \eta, b) \mapsto (\omega, D, \eta, \eta b)$, where $\omega = \nabla^{\perp} \cdot \mathbf{u}$ is the vorticity and $D = \nabla \cdot \mathbf{u}$ is the divergence. If the domain $\Omega$ has no holes, then the vector Laplace equation has a unique solution. This is the Helmholtz theorem for Euclidean space. If the domain does have holes (or islands) then appropriate boundary conditions are required to eliminate harmonic functions in order to guarantee uniqueness of solutions. In either case, provided with appropriate boundary conditions when necessary, one can uniquely split up the vector field $\mathbf{u}$ into potentials via the Helmholtz decomposition

$$
\mathbf{u} = \nabla^{\perp} \psi + \nabla \chi, \quad (13)
$$

where $\psi$ is the stream function and $\chi$ is the velocity potential. Uniqueness of solutions to the Laplace equation is required to reconstruct $\mathbf{u}$, since the potentials satisfy Poisson's equations

$$
\begin{align}
\omega &= \nabla^{\perp} \cdot (\nabla^{\perp} \psi) &&= \Delta \psi, \nonumber \\
D &= \nabla \cdot (\nabla \chi) &&= \Delta \chi. \tag{14}
\end{align}
$$

This permits us to write the velocity field in terms of vorticity and divergence as
follows

$$
\mathbf{u} = \nabla^{\perp} \Delta^{-1} \omega + \nabla \Delta^{-1} D, \quad (15)
$$

which means that changing variables in the Hamiltonian amounts to substitut-
ing (15) for the velocity. The first term in (15) is analogous to the two-dimensional
version of the Biot-Savart law, but since the fluid is compressible it does not
determine the full velocity field. The potential buoyancy variable $r = b/\eta$ is introduced to remove the fraction terms $\pm b_x/\eta$ and $\pm b_y/\eta$ in the Poisson bracket (12). This is convenient from a notation point of view, but also means that all the variables that appear inside the $J$ matrix in the Poisson bracket are “potential variables”, i.e., potential vorticity and potential buoyancy. To change variables in the Cartesian coordinate case, we sandwich the Poisson bracket of Eq. (12) with the functional Jacobian derivative operator and its $L^2$-adjoint as follows

$$ \begin{pmatrix} -\partial_y & \partial_x & 0 & 0 \\ \partial_x & \partial_y & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} 0 & -q & \partial_x & -\frac{b_x}{\eta} \\ q & 0 & \partial_y & -\frac{b_y}{\eta} \\ \partial_x & \partial_y & 0 & 0 \\ \frac{b_x}{\eta} & \frac{b_y}{\eta} & 0 & 0 \end{pmatrix} \begin{pmatrix} \partial_y & -\partial_x & 0 & 0 \\ -\partial_x & -\partial_y & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix} \quad (16) $$

Similar manipulations can be performed for different coordinate systems, provided that one takes an orthonormal coordinate basis. For the sphere this would be a latitude-longitude basis. Performing this change of coordinates in the bracket (12) leads to

$$ \{F, G\} = - \int_{\Omega} \begin{pmatrix} \delta F/\delta\omega \\ \delta F/\delta D \\ \delta F/\delta\eta \\ \delta F/\delta b \end{pmatrix} \nabla \cdot \underbrace{\begin{pmatrix} q \times -q & 0 & r \times \\ q & q \times 1 & r \\ 0 & -1 & 0 \\ r \times -r & 0 & 0 \end{pmatrix}}_{=J(\omega, D, \eta, b)} \nabla \cdot \begin{pmatrix} \delta G/\delta\omega \\ \delta G/\delta D \\ \delta G/\delta\eta \\ \delta G/\delta b \end{pmatrix} d\mu. \quad (17) $$

This bracket formulation has a number of interesting properties. First of all, the differential operators can be factored out of the matrix upon introducing the convention $\times \nabla = \nabla^{\perp}$. Secondly, the matrix only features the potential vorticity $q$ and the potential buoyancy $r$. The kernel of this bracket is the kernel of the matrix operator $J(\omega, D, \eta, b)$, which is skew-symmetric in both $q$ and $r$. The most interesting property in our opinion is that this bracket governs a number of two-dimensional fluid models through submatrices of $J$.

The full four by four setting corresponds to the thermal rotating shallow water equations. If there is no underlying rotation, one simply adapts the definition of the potential vorticity $q$ to obtain the thermal shallow water equations. When there are no buoyancy variations, the bracket can be restricted to the three by three case, which corresponds to the rotating shallow water equations. Again, upon adapting the potential vorticity variable $q$, one can obtain the shallow water equations in case rotation is absent. In the incompressible case where buoyancy variations still play a role, the divergence is zero. Then one can use the submatrix consisting of the (1, 1), (1, 4), (4, 1), (4, 4) elements, which corresponds to the thermal rotating Euler equations. If buoyancy variations do not play a role, the matrix is simply the (1, 1) element, which can describe the two-dimensional rotating Euler equations and the quasi-geostrophic (QG) equations. Important to note is that the thermal
<table><caption>**Table 1** The models that can be described by submatrices of the Poisson bracket (17)</caption><thead><tr><th></th><th>Incompressible</th><th></th><th>Compressible</th></tr></thead><tbody><tr><td>1.</td><td>Euler</td><td>5.</td><td>Shallow water</td></tr><tr><td>2.</td><td>Rotating Euler</td><td>6.</td><td>Rotating shallow water</td></tr><tr><td>3.</td><td>QG</td><td>7.</td><td>Thermal shallow water</td></tr><tr><td>4.</td><td>Thermal rotating Euler</td><td>8.</td><td>Thermal rotating shallow water</td></tr></tbody></table>

QG (TQG) model described in Warneford and Dellar [2013], Zeitlin [2018], Holm et al. [2021] does not fit into the Lie-Poisson bracket formulation because its (1, 1) position features $q-b$ rather than just $q$, see Holm et al. [2021]. This is a result of the fact that the TQG model is derived as a perturbation around thermal geostrophic balance, rather than around geostrophic balance. To summarise, we list all models that can be described by (17) in order of complexity in Table 1.

The transition from incompressible models to compressible models is a steep increase in complexity, since the compressible models involve an additional two equations compared to the incompressible case. The divergence and the depth variable are always paired together, since changes in $D$ imply changes in $\eta$ and vice-versa. This bracket is particularly useful in the explicit computation of Casimirs, which is shown in Sect. 3. It also has other uses. In ocean dynamics, gravity waves propagate at speeds that are orders of magnitude higher than the typical flow velocity. If the Rossby number, the Froude number, the wave amplitude and the stratification parameter satisfy $O(Ro) = O(Fr) = O(\alpha) = O(s)$, then one can derive the thermal geostrophic balance condition. This condition provides an algebraic relation for the balanced velocity field. The balanced velocity field is divergence free. Since the bracket (17) features the divergence variable explicitly, it is natural to perform an asymptotic expansion in a small parameter $\epsilon$ where

$$ 
\begin{aligned}
\omega &= \omega_0 + \epsilon\omega_1 + o(\epsilon), \\
D &= \epsilon D_1 + o(\epsilon).
\end{aligned} 
\qquad (18) 
$$

The expansion (18) applied to the thermal rotating shallow water Lagrangian (1) and truncated at order $o(1)$ yields the thermal extension of the L1 model of Salmon [1983]. This derivation is shown in detail in Holm et al. [2021].

## 3 Explicitly Determining the Casimirs

The bracket (17) is particularly helpful in the explicit computation of the Casimirs since the variables are all independent scalars. We are looking for functionals $C(\omega, D, \eta, b)$ such that $\{F, C\} = 0$ for any functional $F$. The following computations follow a procedure of step-by-step elimination. Expanding the bracket and
requiring {$F, C$} = 0 for any $F$ yields the equation

$$ 
\begin{aligned}
0 = \{F, C\} = & - \int_{\Omega} \frac{\delta F}{\delta \omega} \nabla \cdot \left( q \nabla^{\perp} \frac{\delta C}{\delta \omega} - q \nabla \frac{\delta C}{\delta D} + r \nabla^{\perp} \frac{\delta C}{\delta b} \right) \\
& + \frac{\delta F}{\delta D} \nabla \cdot \left( q \nabla \frac{\delta C}{\delta \omega} + q \nabla^{\perp} \frac{\delta C}{\delta D} + \nabla \frac{\delta C}{\delta \eta} + r \nabla \frac{\delta C}{\delta b} \right) \\
& - \frac{\delta F}{\delta \eta} \nabla \cdot \left( \nabla \frac{\delta C}{\delta D} \right) \\
& + \frac{\delta F}{\delta b} \nabla \cdot \left( r \nabla^{\perp} \frac{\delta C}{\delta \omega} - r \nabla \frac{\delta C}{\delta D} \right) d\mu.
\end{aligned} 
\quad (19) $$

Since $F$ is an arbitrary functional, we can solve for this functional equation per variational derivative of $F$ and get the explicit form of the Casimir s by a process of elimination. In fact, the third line of (19), the one that features $\delta F/\delta \eta$, will be our starting point in the explicit computations. The third line implies that $C$ may be at most linear in $D$ with constant coefficients, since only then the variational derivative of $C$ with respect to $D$ is constant in space. Applying the gradient leads to zero, meaning that the third term vanishes under the assumption of linear dependence of $C$ on mass density $D$. If the variational derivative of $C$ with respect to $D$ is not constant with respect to space, then this term does not vanish, as $\Delta D$ is not necessarily zero. So at this stage we know that

$$ C(\omega, D, \eta, b) = \int_{\Omega} \gamma D + f(\omega, \eta, b) d\mu, \quad (20) $$

where $\gamma \in \mathbb{R}$ is a constant and $f(\omega, \eta, b)$ is a function that is determined next. Since we have established that $C$ must be linear in $D$, all terms that involve variational derivatives of $C$ with respect to $D$ vanish. Simplifying (19) accordingly, we obtain

$$ 
\begin{aligned}
0 = \{F, C\} = & - \int_{\Omega} \frac{\delta F}{\delta \omega} \nabla \cdot \left( q \nabla^{\perp} \frac{\delta C}{\delta \omega} + r \nabla^{\perp} \frac{\delta C}{\delta b} \right) \\
& + \frac{\delta F}{\delta D} \nabla \cdot \left( q \nabla \frac{\delta C}{\delta \omega} + \nabla \frac{\delta C}{\delta \eta} + r \nabla \frac{\delta C}{\delta b} \right) \\
& + \frac{\delta F}{\delta b} \nabla \cdot \left( r \nabla^{\perp} \frac{\delta C}{\delta \omega} \right) d\mu.
\end{aligned} 
\quad (21) $$

The third line of (21), the one that features $\delta F/\delta b$, vanishes if $C$ depends linearly on $\omega$. In this case we do not have to insist on constant coefficients. If the coefficient of $\omega$ is an arbitrary differentiable function $\Psi(r)$ of the potential buoyancy, i.e., $\delta C/\delta \omega = \Psi(r)$, we have

$$ \nabla \cdot \left( r \nabla^{\perp} \frac{\delta C}{\delta \omega} \right) = \nabla \cdot \left( r \nabla^{\perp} \Psi(r) \right) = \Psi'(r) \left( \nabla r \cdot \nabla^{\perp} r \right) + r \left( \nabla \cdot \nabla^{\perp} \right) \Psi(r) = 0, \quad (22) $$
because $\nabla r, \nabla^{\perp} r$ and $\nabla, \nabla^{\perp}$ are orthogonal. Here $\Psi'(r)$ denotes the derivative of $\Psi$ with respect to its argument. Recall that the potential buoyancy is defined as $r = b/\eta$. At this stage, we know that $C$ must have the form

$$C(\omega, D, \eta, b) = \int_{\Omega} \gamma D + \omega \Psi(r) + g(b, \eta) \, d\mu. \qquad (23)$$

The next step is to determine $g(b, \eta)$. Knowing the explicit form of the variational derivative of $C$ with respect to $\omega$, we can simplify (21) further to obtain

$$0 = \{F, C\} = - \int_{\Omega} \frac{\delta F}{\delta \omega} \nabla \cdot \left( q \nabla^{\perp} \Psi(r) + r \nabla^{\perp} \frac{\delta C}{\delta b} \right) + \frac{\delta F}{\delta D} \nabla \cdot \left( q \nabla \Psi(r) + \nabla \frac{\delta C}{\delta \eta} + r \nabla \frac{\delta C}{\delta b} \right) \, d\mu. \qquad (24)$$

We focus on the term that features $\delta F/\delta \omega$. We have an explicitly constructed term and the variational derivative of $C$ with respect to $b$. By the same argument as in the previous step, we know $\nabla \cdot (r \nabla^{\perp}(\delta C/\delta b))$ vanishes if the variational derivative of $C$ with respect to $b$ is a differentiable function of $r$. So, let $\Phi(r) = \Phi(b/\eta)$ be this differentiable function. Since a variational derivative of $\Phi(b/\eta)$ with respect to $b$ produces a factor of $1/\eta$, we introduce the term $\eta\Phi(b/\eta)$ into $C$. So at this stage, after the three steps, we have an expression for $C$ in terms of two arbitrary differentiable functions $\Phi$ and $\Psi$

$$C(\omega, D, \eta, b) = \int_{\Omega} \eta \Phi(r) + \eta q \Psi(r) + \gamma D \, d\mu \qquad (25)$$

where $\gamma \in \mathbb{R}$ is a constant and $\eta q = \omega + \frac{1}{R_0} f$. Going back to (24), we can verify whether $C(\omega, D, \eta, b)$ is indeed the family of Casimir s of the Poisson bracket (17). This means that the terms multiplying variational derivatives of $F$ must vanish. So we substitute (25) into (24) and compute the term multiplied by $\delta F/\delta \omega$

$$\begin{align} & \nabla \cdot \left( q \nabla^{\perp} \frac{\delta C}{\delta \omega} - q \nabla \frac{\delta C}{\delta D} + r \nabla^{\perp} \frac{\delta C}{\delta b} \right) \nonumber \\ &= \nabla \cdot \left( q \nabla^{\perp} \Psi(r) + r \nabla^{\perp} (\Phi'(r) + q \Psi'(r)) \right) \nonumber \\ &= \Phi''(r)(\nabla r \cdot \nabla^{\perp} r) + r(\nabla \cdot \nabla^{\perp})\Phi'(r) \tag{26} \\ & \qquad + \Psi'(r)(\nabla q \cdot \nabla^{\perp} r) + \Psi'(r)(\nabla r \cdot \nabla^{\perp} q) \nonumber \\ & \qquad + r\Psi''(r)(\nabla r \cdot \nabla^{\perp} q) + r\Psi''(r)(\nabla q \cdot \nabla^{\perp} r) \nonumber \\ &= 0. \nonumber \end{align}$$
In this computation we have used orthogonality of $\nabla$, $\nabla^\perp$ and $\nabla r$, $\nabla^\perp r$ and skew-symmetry, i.e., $\nabla r \cdot \nabla^\perp q = -\nabla q \cdot \nabla^\perp r (= -r_x q_y + q_x r_y)$ for Cartesian coordinates). It remains to check whether the term multiplied by $\delta F/\delta D$ also vanishes. In the following computation, we suppress the dependence of $\Phi$ and $\Psi$ on $r$ for notational convenience. We compute from the second line in (19)

$$ 
\begin{aligned}
& \nabla \cdot \left( q \nabla \frac{\delta C}{\delta \omega} + q \nabla^\perp \frac{\delta C}{\delta D} + \nabla \frac{\delta C}{\delta \eta} + r \nabla \frac{\delta C}{\delta b} \right) \\
&= \nabla \cdot \left( q \nabla \Psi + \nabla (\Phi - r \Phi' - r q \Psi') + r \nabla (\Phi' + q \Psi') \right) \quad (27) \\
&= \nabla \cdot \left( q \nabla \Psi - r \nabla (q \Psi') - q \Psi' \nabla r + r \nabla (q \Psi') \right) \\
&= 0.
\end{aligned}
$$

The computation is a sequence of applying the identity $\nabla f(r) = f'(r) \nabla r$. In the first step we have applied this to all the terms that involve $\Phi$. Performing the same manipulations on $\Psi$ subsequently, implies the result. Note that in (24) we required cancellations of terms orthogonal to the terms that are required to cancel in (27). The reason that the perpendicular gradient of the variational derivative of $C$ with respect to $\eta$ does not appear in (24) is because it is trivially zero due to orthogonality of $\nabla$ and $\nabla^\perp$. Hence we can conclude that (25) is the complete description of the Casimirs for the bracket (17) and thus also for (11). We can repeat the computation for each of the models described in Table 1 to obtain the corresponding Casimirs. This is summarised in Table 2. Here one can see that in presence of thermal effects enstrophy is no longer a Casimir.

<table><caption>**Table 2** The Casimirs of the models that can be described by submatrices of the Poisson bracket (17)</caption><thead><tr><th></th><th>Incompressible</th><th>Casimirs</th></tr></thead><tbody><tr><td>1.</td><td>Euler</td><td>C = ∫<sub>Ω</sub> Φ(q) dμ</td></tr><tr><td>2.</td><td>Rotating Euler</td><td>C = ∫<sub>Ω</sub> Φ(q) dμ</td></tr><tr><td>3.</td><td>QG</td><td>C = ∫<sub>Ω</sub> Φ(q) dμ</td></tr><tr><td>4.</td><td>Thermal rotating Euler</td><td>C = ∫<sub>Ω</sub> Φ(r) + qΨ(r) dμ</td></tr><tr><td></td><td>Compressible</td><td>Casimirs</td></tr><tr><td>5.</td><td>Shallow water</td><td>C = ∫<sub>Ω</sub> ηΦ(q) + γD dμ</td></tr><tr><td>6.</td><td>Rotating shallow water</td><td>C = ∫<sub>Ω</sub> ηΦ(q) + γD dμ</td></tr><tr><td>7.</td><td>Thermal shallow water</td><td>C = ∫<sub>Ω</sub> ηΦ(r) + ηqΨ(r) + γD dμ</td></tr><tr><td>8.</td><td>Thermal rotating shallow water</td><td>C = ∫<sub>Ω</sub> ηΦ(r) + ηqΨ(r) + γD dμ</td></tr></tbody></table>
# 4 Conclusion

In geometric approaches to fluid dynamics one often may exploit constructive methods to infer conservation laws via Noether’s theorem Abraham and Marsden [1978], Marsden and Weinstein [1983], Marsden and Ratiu [2013], Holm et al. [2009]. Casimirs identify conservation laws that arise as the kernel of the Poisson bracket. We provided an explicit method of determining these Casimir functionals for two-dimensional fluid dynamics by means of two changes of variables for the thermal rotating shallow water equations. We formulated the equations in vector invariant form by using an important vector calculus identity. The Poisson bracket corresponding to this vector invariant form is convenient for further manipulation. By changing coordinates from velocity to vorticity and divergence, we derive a Poisson bracket that only involves (skew) gradients and divergences. By means of this formulation, it is a systematic computation to obtain the Casimirs. The computations were performed for arbitrary coordinate systems, which means that the above computations can easily used for domains such as the sphere.

**Acknowledgments** The authors are grateful for many fruitful discussions with Paolo Cifani, Sagy Ephrati, Arnout Franken and Darryl Holm. We thank the anonymous reviewer for their valuable comments.

# References

Ralph Abraham and Jerrold Marsden. *Foundations of mechanics*, volume 36. Benjamin/Cummings Publishing Company Reading, Massachusetts, 1978.

Colin Cotter and Darryl Holm. A variational formulation of vertical slice models. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 469(2155):20120678, 2013.

Darryl Holm and Bruce Long. Lyapunov stability of ideal stratified fluid equilibria in hydrostatic balance. *Nonlinearity*, 2(1):23, 1989.

Darryl Holm and Erwin Luesink. Stochastic wave-current interaction in thermal shallow water dynamics. *Journal of Nonlinear Science*, 31(2):1-56, 2021.

Darryl Holm, Jerrold Marsden, and Tudor Ratiu. The Euler-Poincaré equations and semidirect products with applications to continuum theories. *Advances in Mathematics*, 137(1):1-81, 1998.

Darryl Holm, Tanya Schmah, and Cristina Stoica. *Geometric mechanics and symmetry: from finite to infinite dimensions*, volume 12. Oxford University Press, 2009.

Darryl Holm, Erwin Luesink, and Wei Pan. Stochastic mesoscale circulation dynamics in the thermal ocean. *Physics of Fluids*, 33(4):046603, 2021.

Robert Kraichnan. Inertial ranges in two-dimensional turbulence. *The Physics of Fluids*, 10(7): 1417-1423, 1967.

Erwin Luesink. Stochastic geometric mechanics of thermal ocean dynamics. *PhD thesis*, 2021.

Jerrold Marsden and Tudor Ratiu. *Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems*, volume 17. Springer Science & Business Media, 2013.

Jerrold Marsden and Alan Weinstein. Reduction of symplectic manifolds with symmetry. *Reports on mathematical physics*, 5(1):121-130, 1974.
Jerrold Marsden and Alan Weinstein. Coadjoint orbits, vortices, and clebsch variables for incompressible fluids. *Physica D: Nonlinear Phenomena*, 7(1-3):305-323, 1983.

Rick Salmon. Practical use of Hamilton's principle. *Journal of Fluid Mechanics*, 132:431-444, 1983.

Emma Warneford and Paul Dellar. The quasi-geostrophic theory of the thermal shallow water equations. *Journal of Fluid Mechanics*, 723:374-403, 2013.

Vladimir Zeitlin. Geophysical fluid dynamics: understanding (almost) everything with rotating shallow water models. Oxford University Press, 2018.

**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_85.png" /></figure>
# Correlated Structures in a Balanced Motion Interacting with an Internal Wave

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

Check for
updates

Igor Maingonnat, Gilles Tissot, and Noé Lahaye

**Abstract** Characterizing the loss of coherence of an internal tide propagating through mesoscale turbulence has been a major challenge in oceanography, particularly due to its implications for the interpretation of satellite data. In this paper, we intend to study the correlations between a balanced motion and the incoherent part of a wave in an idealised configuration. We introduce a new modal decomposition technique, named broad-band proper orthogonal decomposition (BBPOD), which consists in performing a proper orthogonal decomposition (POD) on complex demodulated variables. After connecting BBPOD to the standard SPOD, we show that BBPOD, coupled with the extended POD technique enables us to associate the principal components of the incoherent field to the slow flow structures responsible of this loss of coherence through triadic interactions with the incident wave.

**Keywords** Internal tide interactions · Spectral proper orthogonal decomposition · Broadband proper orthogonal decomposition

## 1 Introduction

Internal tides, generated by interactions between the barotropic tide and topographic features such as ridges or continental slopes, are ubiquitous in the ocean, playing a crucial role in vertical mixing and energy transport. They propagate over large distances, encountering regions with energetic mesoscale turbulence, and they lose their fixed phase relationship with the astronomical forcing, a phenomenon known as incoherence. This loss of coherence, highly unpredictable, complicates for example our ability to disentangle internal tides and low-frequency turbulent signals from satellite data (Richman et al. 2012).

I. Maingonnat (✉) · G. Tissot · N. Lahaye

INRIA Rennes Bretagne Atlantique, IRMAR – UMR CNRS 6625, Rennes, France

e-mail: igor.maingonnat@inria.fr
To face these difficulties, there have been various studies to better understand and predict the impact of a balanced (turbulent) jet on the inertia-gravity wave field properties. For instance, Savva and Vanneste (2018) and Ward and Dewar (2010) studied the internal tide scattering. Ponte and Klein (2015) examined the incoherence time scales for different turbulent configurations, and Dunphy et al. (2017) quantified the interaction terms via vertical-mode projection of the linearized Boussinesq equations.

We consider in the present paper a data-driven approach to study from idealised numerical simulations, the structures of the jet, which are correlated with the incoherent contributions of the wave field. Extracting coherent structures from a data set can be performed in the spectral domain by spectral proper orthogonal decomposition (SPOD) (Towne et al. 2018). Some attempts to connect these reduced features with non-linear interactions have been proposed for instance in Karban et al. (2022) by identifying the contribution of the non-linear term correlated with the SPOD mode through extended proper orthogonal decomposition (EPOD, Boree (2003)). Unfortunately, non-linear forcing is a quantity more difficult to interpret and associate with physical mechanisms than a pressure or velocity field for instance. We propose instead a new *broad-band proper orthogonal decomposition (BBPOD)*, derived from the SPOD. By taking advantage of the strong time scale separation between the two dynamics and considering small-amplitude wave, this formulation allows us to connect the non-linear interactions between the slow motions of the jet and the incoherent contributions of the wave through the EPOD method. The non-linear interactions are here understood as triadic interactions with the slow motion that will lead to a broadband frequency structure of the incoherent wave field.

The plan is as follows. We will begin by describing the model used for our simulations in Sect. 2. In Sect. 3.1, the original SPOD method is reviewed and a connection to the proposed BBPOD method is made (Sect. 3.2). Sections 4 and 5 will summarize the study and bring some perspectives.

## 2 Model

The propagation of internal tides through a nearly-balanced jet is examined in a one layer rotating shallow water (RSW) model. The equations are non-dimensionalized as follows. The characteristic time-scale is the inertial time $T = f^{-1}$, inverse of the Coriolis frequency. The reference length scale $L$ is chosen to be of the order of the jet thickness, taken equal to the first Rossby radius of deformation $R_d$, such that the Burger number $B_u = R_d^2/L^2$ is equal to one. A beta-plane approximation, with parameter $\beta$, accounts for the effect of rotation, and a radiative damping term $\alpha$ is added on the continuity equation, (following e.g. Brunet and Vautard 1996). Adequat artificial hyperviscosity, consistent with energy dissipation and conservation of angular momentum, is also used (see Ochoa et al. 2011).
Let $\Omega$ be a bounded subset of $\mathbb{R}^2$, the equations without viscosity, defined on $\Omega \times \mathbb{R}_+$, are

$$ \partial_t h + B_u \mathrm{div} \mathbf{v} + R_o((\mathbf{v} \cdot \nabla)h + h \mathrm{div} \mathbf{v}) = -\alpha h \quad (1) $$

$$ \partial_t \mathbf{v} + (1 + \beta y) \mathbf{v}^{\perp} + R_o(\mathbf{v} \cdot \nabla) \mathbf{v} = -\nabla h, \quad (2) $$

where $\mathbf{v} = (u, v)$ is the horizontal velocity, $\mathbf{v}^{\perp} = (-v, u)^T$, $h$ is the sea surface height (SSH) anomaly and $R_o = U/(fL)$ is the Rossby number.

A numerical simulation of system (1) and (2) of a plane wave interacting with a zonal jet has been performed. Equations have been discretized by a pseudo-spectral method and with a Runge-Kutta time-scheme using the code Dedalus (Burns et al. 2020). The domain $\Omega$ is a doubly periodic rectangular domain of size $[0, 20] \times [-20, 20]$ discretized on a $128 \times 512$ grid. The simulation is initialized with an eastward zonal jet at geostrophic equilibrium with a small perturbation to trigger the instabilities. During the whole experiment, a northward propagating plane wave with frequency fixed at $\omega = 2$ is generated in a nudging layer in the South of the domain, and an eastward wind forcing is applied to maintain the balanced current. A spin-up phase allows the jet to evolve toward a statistically stationary state. The experiment continues with snapshots stored every 0.1 wave period, such that a sufficiently long series $q = (u, v, h)^T$ of 12,000 snapshots are collected, representing nearly 4000 wave periods. Note that a sufficient sampling in time is required for extracting the wave field by filtering in Sect. 3.2. An example of snapshot associated with one run of the simulation is shown in Fig. 1.

<figure><img src="image_87.png" /><figcaption>Fig. 1 A snapshot taken from the RSW simulation. On the left the vorticity field and on the right the internal wave SSH anomaly $h_\omega$, extracted by bandpass filtering centered around the wave frequency $\omega = 2$</figcaption></figure>
This simple model can sustain waves interacting with a turbulent flow. Northward rotating shallow water waves are inertia gravity waves satisfying the dispersion relation

$$k_y^2 = \frac{\omega^2 - (1 + \beta y)^2}{B_u}, \qquad (3)$$

with $k_y$ denoting the spatial wavenumber in the meridional direction, which is function of the tidal frequency $\omega$.

# 3 Methods

The methods to extract correlated structures in a flow field are described in this section. In the following we denote our state vector by $q = (u, v, h)^T$.

## 3.1 Spectral Proper Orthogonal Decomposition

Spectral proper orthogonal decomposition is an extension of POD methods (see Lumley 1967), and aims at extracting coherent structures in spectral space from numerical or experimental data. The data are assumed to be issued from a statistically stationary random process, which is verified thanks to the hypothesis of small-amplitude waves, such that the balanced flow is marginally impacted by the wave propagation. The modes obtained from this decomposition are space-time uncorrelated from each other.

We interpret our state vector as a zero-mean second order random process $q$ (which can be obtained by subtracting the time averaged field beforehand) indexed over $\Omega \times \mathbb{R}^+$, a subset of spatio-temporal variables. Consider $\hat{q}(\omega)$ the Fourier transform of the process at angular frequency $\omega$. We assume that each realization of $\hat{q}(\omega)$ belongs to $L_W^2(\Omega, \mathbb{C}^3) = \{g : \Omega \to \mathbb{C}^3, \int_{\Omega} g^*(z)W(z)g(z) dz < \infty\}$, the ensemble of square integrable functions in space (periodic at boundary) relatively to a positive definite weight matrix $W: \Omega \to M_n(\mathbb{R})$. The superscript $*$ denotes the transpose-conjugate operation. The matrix $W$ is chosen such that the $L^2$-norm approximates the energy of the RSW model (1) and (2):

$$\|q\|_{\mathbf{W}}^2 = \frac{1}{2} \int_{\Omega} (1 + \lambda \bar{h})(u^2 + v^2) dxdy + \frac{1}{2B_u} \int_{\Omega} h^2 dxdy, \qquad (4)$$

where $\lambda$ refers to the deviation of the isopycnal taken equal to $R_o/B_u$ and $\bar{h}$ corresponds to the temporal-mean of SSH anomaly. The term $1 + \lambda \bar{h}$ is assumed to remain strictly positive (which is necessary for ensuring positiveness of the norm), since we do not consider zero or negative sea-surface height.
Let us recall briefly the basic principle of SPOD. By stationarity the auto-correlation function $C(x, y, t, x', y', t') = E(q(x, y, t) \otimes q^*(x', y', t'))$ satisfies:

$$C(x, y, t, x, y, t') = C(x, y, x', y', t-t'), \quad (5)$$

where $\otimes$ denotes the dyadic product of two vectors in $\mathbb{R}^3$. The objective of SPOD method is therefore to find a deterministic function $\psi$ satisfying the Fredholm equation:

$$\int_{\Omega} \mathbf{S}(x, y, x', y', \omega) \mathbf{W}(x', y') \psi(x', y', \omega) \, dx'dy' = \lambda(\omega) \psi(x, y, \omega), \quad (6)$$

where $\mathbf{S}(x, y, x', y', \omega) = \int_{\mathbb{R}} \mathbf{C}(x, y, x', y', \tau) e^{-i\omega\tau} d\tau$ is the cross spectral density matrix (CSD).

The Karhunen-Loeve theorem (Loève 1955) states that $\forall \omega \in \mathbb{R}$, Eq. (6) has an infinite number of solutions $(\lambda_j(\omega), \psi_j(\omega))_{j=1}^{\infty}$ such that: $(\psi_j(\omega))_{j=1}^{\infty}$ is an orthonormal basis in $L^2_{\mathbf{W}}(\Omega, \mathbb{C}^3)$ where we can expand the Fourier transform of the field into structures uncorrelated from each other:

$$\hat{\mathbf{q}}(x, y, \omega) = \sum_{j=1}^{\infty} a_j(\omega) \psi_j(x, y, \omega) \quad (7)$$

with $E(a_j(\omega)a_{j'}(\omega')) = \lambda_j(\omega)\delta_{j,j'}(\omega - \omega')$ and $\forall j, \lambda_j(\omega) \ge 0$. Moreover, a truncated expansion at order $n$ will maximize the mean energy (defined by the norm $\|\cdot\|^2_{\mathbf{W}}$), compared to any other decomposition of the same order, and SPOD provides an optimal decomposition of the CSD:

$$\mathbb{E}(\|\hat{\mathbf{q}}(\omega)\|^2_{\mathbf{W}}) = \sum_{j=1}^{\infty} \lambda_j(\omega) \quad (8)$$

$$\mathbf{S}(x, y, x', y', \omega) = \sum_{j=1}^{\infty} \lambda_j(\omega) \psi_j(x, y, \omega) \psi_j^*(x', y', \omega). \quad (9)$$

Therefore, by solving the Fredholm equation at the tidal frequency, we are able to separate the fast from the slow component of the process by a spectral decomposition, expanded onto an orthonormal basis optimal energetically.

It can be remarked that when there is an homogeneous direction, it is possible to compute the Fredholm equation in spectral space with respect to this direction. Each wave-number can be computed independently, and each corresponding Fourier component is solution of the original Fredholm problem. For more complete description of SPOD (see Towne et al. 2018; Schmidt and Colonius 2020). In our case, the domain is homogeneous in the $x$ direction, but we still consider the 2D
problem for the physical relevance of the extended BBPOD problem presented in Sect. 3.2.

## 3.2 Broadband Proper Orthogonal Decomposition

The BBPOD consists in estimating by complex-demodulation (see Godfrey 1965) the correlation tensor of a band-passed-filtered signal, and its eigenfunctions. The associated algorithm is presented in this section. It will be compared with the Welch's method used in SPOD (see Welch 1967; Towne et al. 2018) and it will be shown that under some hypotheses both algorithms are equivalent. It can be noted that in Welch (1967), the connection with complex demodulation is briefly mentioned, and it is leveraged here with the computation of the eigenvectors and eigenvalues of the correlation tensor for a POD decomposition. We highlight by this equivalence, that due to windowing in the Welch's method, the features extracted by SPOD possess a spectral component with a thick frequency band. For the Broadband POD algorithm, this frequency band is explicit, and chosen through the definition of a filter. In the RSW model, the incoherent wave field possesses a broad band structure and Broadband POD allows us to obtain a complete decomposition of this field.

### 3.2.1 Complex Demodulation of the Wave Field

Given a temporal series $x_t$, the principle of complex demodulation consists in the following computation:

$$x_d = \mathcal{L}(x_t e^{-i\omega t}), \qquad (10)$$

where $x_d$ is called the complex demodulated signal, $\mathcal{L}$ is a low-pass filter and $\omega$ is the demodulated frequency. Complex demodulation enables to extract the slowly varying amplitude envelop and phase deviations of the wave field. These variations are associated with wave incoherence.

First, we decompose the process $q(x, y, t)$ in a sum of $q_j(x, y, t)$ and $q_\omega(x, y, t)$, representing the jet and wave contributions, respectively. The wave field can be expressed as:

$$q_{\omega}(x, y, t) = q_{d}(x, y, t)e^{i\omega t}, \qquad (11)$$

where $\omega$ is the tidal frequency and $q_d$ is a complex field slowly varying in time, accounting for the incoherence. The latter is then extracted by the complex demodulation of $q$ (10), with a filter designed such that it isolates the slow variations of the background flow. As a consequence, through the $e^{-i\omega t}$ shift, the spectral-band of the extracted wave contains all effects of triadic interactions between the
slow motions and the coherent wave. The jet contribution $q_j$ is obtained by directly filtering $q$ (zero frequency), including here the time-average. In this study, we assume that $q$ is dominated by the superimposition of $q_j$ and $q_\omega$, which interact non-linearly. Note that for simplicity, a uniform time sampling of the data was chosen but this can be adapted to a non-uniform one, e.g. by interpolation on a regular grid or using filtering algorithms to irregular sampling.

### 3.2.2 Link with SPOD

For making the link with SPOD to estimate the CSD of a signal $x_t$, we write the temporal filter in a discrete form, with a time spacing $\Delta t$ and $t_k = k\Delta t$. A wide class of linear filter can be expressed as the convolution:

$$ \mathcal{L}(x_t) = \sum_{i=-m}^{m} b_i x_{t-i}, \qquad (12) $$

where $(b_i)_{-m \le i \le m}$ are the discrete coefficients of the filter. Then,

$$ (\mathcal{L}(x_t e^{-i\omega t}))_j = \sum_{k=-m}^{m} b_k x_{j-k} e^{-i\omega t_{j-k}} \qquad (13) $$

$$ = \sum_{k=j-m}^{j+m} b_{j-k} x_k e^{-i\omega t_k}. \qquad (14) $$

In the Welch method, $x_t$ is subdivided into possibly overlapping blocks of size N and with an overlap $N_o$. A Fast Fourier Transform is performed on each windowed block to extract the Fourier component at the tidal frequency, denoted $X_\omega^l$ where $l$ is the block index. We define

$$ X_{\omega}^{l} = \sum_{k=0}^{N} x_{k+l(N-N_o)} W_k e^{-i\omega t_k}, \qquad (15) $$

where $W_k$ is a window function. By changing variable $k' = k + l(N - N_o)$, we obtain

$$ X_{\omega}^{l} = e^{i\omega t_{l}(N-N_{o})} \sum_{k'=l(N-N_{o})}^{N+l(N-N_{o})} x_{k'} W_{k'-l(N-N_{o})} e^{-i\omega t_{k'}}. \qquad (16) $$

Assuming that the window function is symmetrical in the middle of each block (which is verified for most windows used in the literature), i.e. $W_k = W_{N-k}$, Eq. (16) gives
$$X_{\omega}^{l} = e^{i\omega t_{l(N-N_{o})}} \sum_{k'=l(N-N_{o})}^{N+l(N-N_{o})} W_{N+l(N-N_{o})-k'x_{k'}} e^{-i\omega t_{k'}}. \quad (17)$$

Finally, by choosing $W_k = b_k$ which sets $m = \frac{N}{2}$, relation (14) yields

$$X_{\omega}^{l} = e^{i\omega t_{l(N-N_{o})}} (\mathcal{L}(x_{t}e^{-i\omega t}))_{\frac{N}{2}+l(N-N_{o})}. \quad (18)$$

Consequently, up to a phase, the Welch method corresponds to the computation of the complex demodulated signal sampled every $N - N_o$. The phase shift cancels when computing the correlation over $N_b$ blocks:

$$\sum_{l=0}^{N_b} X_{\omega}^{l} X_{\omega}^{l*} = \sum_{l=0}^{N_b} (\mathcal{L}(x_t e^{-i\omega t}))_{\frac{N}{2}+l(N-N_o)} (\mathcal{L}^*(x_t e^{-i\omega t}))_{\frac{N}{2}+l(N-N_o)}. \quad (19)$$

The CSD can thus be obtained by computing the correlation over $N_b$ snapshots sampled every $N - N_o$ of the complex demodulated signal, which is proving the equivalence between the two methods for appropriate numerical parameters.

As said before, it can be remarked that the window function acts as a filter in the Welch procedure, but without giving an explicit expression of the frequency band. Moreover, if we aim at studying the whole broadband field, the classical SPOD (with its single frequency interpretation when the window function is designed to only estimate a pure harmonic) requires the computation of a set of spatial modes at each discrete frequency in the peak. In comparison, BBPOD computes at only one frequency the dominant modes of the wave field. This has the drawback for the SPOD algorithm to include at each frequency a few modes of nearby frequency due to the convolution with a finite length window, leading to misinterpretation or counting the same mode several times in the reconstruction of the whole wave field. More precisely, BBPOD gathers the first SPOD mode of each frequency providing a simpler representation of the whole wave field and avoids spurious modes in the reconstruction. Another important remark is that the connection of our method with SPOD lies on the fact that the spectral peak is sufficiently narrow corresponding of a sufficiently large window in the Welch method guaranteeing a good estimation of the CSD.

The final step is to compute the eigenvalues and eigenvectors of the CSD, estimated by complex demodulation. The complex demodulated signal at frequency $\omega$ is then decomposed on this BBPOD basis:

$$\mathbf{q}_d(x, y, t) = \sum_j a_j(\omega, t) \psi_j(x, y, \omega), \quad (20)$$

with $a_j(\omega, t) = \int_{\Omega} \mathbf{q}_d^*(x, y, t) \mathbf{W}(x, y) \psi_j(x, y, \omega) dx dy$ slowly varying coefficients.
### 3.2.3 Extended Broadband Proper Orthogonal Decomposition

Non-linear interactions between slow variations of the background flow and the wave induce incoherent wave contributions through triadic interactions. A major interest of BBPOD is that it allows us to study the correlations between the slow perturbations of the jet, and the incoherence of the fast wave field, by extracting the slowly varying complex demodulated amplitudes. To that end, we propose to apply the concept of EPOD in the framework of BBPOD. EPOD, originally presented in Boree (2003), enables to identify the part of a target field correlated with a given POD mode. The target field we first consider is the balanced motion obtained by low-pass filtering. The p-th extended POD mode of the slow motion $q^j$ correlated to the p-th BBPOD coefficient of the wave $a_p^\omega$ is defined by:

$$ \chi_p^j = \frac{\mathbb{E}(q^j a_p(\omega))}{\lambda_p(\omega)}, \qquad (21) $$

where the expectation operator is a temporal average over snapshots. The mode $\chi_p^j$ will be referred to as direct EPOD, and $a_p(\omega)\chi_p^j$ represents the part of the jet correlated with the p-th broadband POD coefficient of the wave.

Complementary to direct EPOD, we define as well an inverse EPOD, applying BBPOD for $\omega = 0$, thus obtaining an orthonormal basis representative of the jet variability, and considering the complex demodulated of the wave contribution $q_d$ as the target field. In this situation, we identify the contribution of the incoherent wave field correlated with the BBPOD coefficients of the jet noted $a_p(\omega = 0)$. Therefore, in the following direct EPOD refer to jet modes while inverse EPOD refer to wave modes.

## 4 Results

This section details the numerical results carried out in this work. The goal is to understand the non-linear interactions between the balanced motion and the wave field by means of the methods presented in Sect. 3, that will extract the correlations between the two dynamics.

Figure 2 shows the real part of the first three energy-scaled BBPOD modes at the tidal frequency, sorted by decreasing energy (4). Specifically here, we plot the scaled quantities $\sqrt{\lambda_j}\psi_j$ to highlight their respective energy contribution, but the analysis was performed with the normalized modes. To define the low pass filter $\mathcal{L}$, a fourth order Butterworth filter has been taken with a frequency cut-off equals to the typical frequencies of the jet, such that the whole wave-field scattered by the jet is captured. The first mode containing the most energy corresponds to the coherent wave field, that is defined here as the part of the wave correlated to the tidal forcing and thus phase-locked in time, since it corresponds to the time-average
<figure><img src="image_88.png" /><figcaption>**Fig. 2** Real part of the first three energy-scaled BBPOD modes $\sqrt{\lambda_j \psi_j}$ associated with $u$ at the tidal frequency. Sponge regions are not shown</figcaption></figure>

<table><caption>**Table 1** Normalized energy of the first 5 modes shown in percentage. Computed in $\Omega$ and in $\Omega^N$. Each eigenvalue is divided by the total variance in $\Omega$ at the tidal frequency $\sum_j \lambda_j^\Omega(f)$</caption><thead><tr><th>Mode</th><th>Ω</th><th>Ω<sup>N</sup></th></tr></thead><tbody><tr><td>1</td><td>57.35</td><td>29.04</td></tr><tr><td>2</td><td>13.81</td><td>12.98</td></tr><tr><td>3</td><td>9.86</td><td>9.51</td></tr><tr><td>4</td><td>3.86</td><td>3.77</td></tr><tr><td>5</td><td>2.34</td><td>2.13</td></tr><tr><td>Sum of incoherence</td><td>42.64</td><td>37.91</td></tr></tbody></table>

of the complex-demodulated signal. Other modes account then for the incoherent structures as space-time decorrelated of the mean. They show deflections in opposite direction, with mode-number $m_x = \pm 2$, and with meridional mode-number $m_y$ expected to be determined by the dispersion relation (3) and the structure of the jet (see Bühler 2014). This is a consequence of the homogeneity in the zonal direction, guaranteed by our idealised set-up, which makes SPOD modes equal to Fourier modes in $x$.

Table 1 shows the normalised energy contained in the first five BBPOD modes and in the incoherence field. It is computed for the domain without sponge regions noted $\Omega = [0, 20] \times [-12, 12]$ and for the subset $\Omega^N = [0, 20] \times [-2, 12]$ representing the region of incoherences.

In the total domain, the coherent plane wave accounts for 57.35% of the energy and the first three modes 81% of the energy. This shows that linear effects dominate the wave propagation, because we are in a configuration for which incoherence is relatively weak. As explained in Ward and Dewar (2010), a more energetic jet interacting with a tide with higher frequency would result in a more energetic incoherent wave field for a RSW simulation. As expected, incoherent modes have nearly all their energy in the upper part while it is rather equally distributed for the coherent mode. In this region, incoherences are dominant due to non-linearities
<figure><img src="image_89.png" /><figcaption>Fig. 3 Cumulative energy $\sum_{j=1}^{N} \lambda_j(\omega) / \sum_{j=1}^{+\infty} \lambda_j(\omega)$ [%] contained in the reconstruction of the incoherent field by BBPOD modes</figcaption></figure>

in the center of the domain, which leads to an increase of the amplitude of the modes with y, even though the coherent mode still contains the most energy. Figure 3 shows the energy contained in the reconstructed incoherent wave field in terms of cumulative energy. The first two incoherent modes account for 58% of the total incoherent energy. For 6 modes, 80% of the energy is recovered in the reconstruction. The cumulative energy can be understood as a normalized RMSE measuring the accuracy of the basis to reconstruct the true incoherent wave field $q_{inc}$

according to $\frac{\mathbb{E}(\|q_{inc} - \sum_{j=1}^{N} (q_{inc} \cdot \psi_j) \psi_j\|_{L^2(\Omega)}^2)}{\mathbb{E}(\|q_{inc}\|_{L^2(\Omega)}^2)} = 1 - \sum_{j=1}^{N} \lambda_j(\omega) / \sum_{j=1}^{+\infty} \lambda_j(\omega)$. It can be remarked that spatial coherence can be computed with the coherence function $\gamma$ through the CSD.

A SPOD decomposition was also performed in Egbert and Erofeeva (2021) for a realistic HYCOM simulation. Similar results were found for the energy captured by the SPOD modes at the M2 frequency, which is similar to our forcing frequency. This encourages to extend this analysis to a more realistic simulation.

Figure 4 depicts the direct EPOD modes of the balanced motion, correlated with the BBPOD modes of the wave. The modes are weighted as explained before. The first extended mode is showing that the part of the jet correlated with the coherent field is the stationary mean flow, which connects the first order statistic of both dynamics. This results is the expression of the linear propagation of the incident wave through the mean flow (in time and in the x-direction). The second and third EPOD modes represent the meanders of the jet, consisting of a vortex train. The modes are defined up to a complex constant. We can see that $\chi_2^j$ is approximately proportional to $\chi_3^j$, suggesting that meandering of the jet generates eastward and westward propagating incoherent perturbations.
<figure><img src="image_90.png" /><figcaption>**Fig. 4** On the left column are presented the real part of the first three BBPOD modes of $u$ at the tidal frequency. On the last column are the corresponding direct EPOD modes of $u$. Modes are scaled as follows: $\sqrt{\lambda_j} \psi_j$</figcaption></figure>

The opposite procedure is done next by inverse EPOD. Figure 5 shows the leading three jet (energy-scaled) BBPOD modes associated with $\omega = 0$ and the associated EPOD modes of the wave. As in direct EPOD, the first mode is showing the correlation between the mean of both dynamics. For the second and third modes,
<figure><img src="image_91.png" /><figcaption>Fig. 5 On the first column are presented the first three BBPOD modes of $u$ of the jet. On the last column are the associated inverse EPOD modes of $u$. Modes are scaled as follows: $\sqrt{\lambda_j \psi_j}$</figcaption></figure>

BBPOD shows jet meandering structures identical to the direct EPOD modes (up to a phase). However, their associated inverse EPOD indicate a standing wave corresponding to the sum of the two direct broadband POD modes that represent left and right deflections, expressing the fact that the EPOD and BBPOD association is not bijective. A given meander will give rise to a wave deflected eastward nearly as much as a wave deflected westward, and their relative contribution is not possible
to predict, but the present result suggests that an approximately equal repartition
would constitute a fair estimate.

Globally, we can infer that the most energetic triadic interactions between the jet and the plane wave are caused by meandering structures of the jet generating an incoherent field of the form of a standing wave in the zonal direction x and propagating northward. The rising of a standing wave suggests in average an equal distribution of the eastward and westward propagating modes. This gives nice perspectives of incorporating this methodology in the context of data assimilation to estimate the incoherent wave field with some knowledge of the slow balanced flow through triadic interactions.

**5 Summary and Perspectives**

This study investigates the correlated processes at stake within the non-linear
interactions of a plane wave propagating through a meso-scale oceanic current. For
that purpose, the method of broadband POD is introduced and connected to extended
POD. A first methodological point is done to show that BBPOD is equivalent to
the common Welch method to estimate the cross spectral density matrix, to which
SPOD are the eigenfunctions. More importantly for our study, connecting these
two methods allows us to understand that the most energetic features of the wave
are obtained by triadic interactions with the most energetic features of the jet. In
our idealised configuration, the wave is deflected in both directions by a given
meander. The spatial periodicity of the meander sets the zonal deflections, giving
us clearer insight on the BBPOD of the wave knowing only the extended modes
of the jet. This method can be seen as complementary to methods based on the
bispectrum to detect triadic interactions in the data. Here we have an a priori on
the nature of the interactions, between a low frequency jet and a fast wave field,
allowing to study directly their correlations and to propose a methodology well
suited for reconstruction. In a future work, we intend to examine the sensitivity of
the results to other jet or wave configurations, and even to richer models. Studying
the implications of these methods to reconstruct by correlation the internal wave
field by the knowledge of the instantaneous slow motion is also envisaged for data
assimilation problematics, either by direct estimation, or for learning an observation
operator associated with a Galerkin model defined in the BBPOD space.

**References**

Boree, Jacques (Aug. 2003). "Extended proper orthogonal decomposition: A tool to analyse correlated events in turbulent flows". In: *Experiments in Fluids* 35, pp. 188–192. DOI: https://doi.org/10.1007/s00348-003-0656-3.
Brunet, Gilbert and Robert Vautard (1996). "Empirical Normal Modes versus Empirical Orthogonal Functions for Statistical Prediction". In: *Journal of Atmospheric Sciences* 53.23, pp. 3468-3489. DOI: https://doi.org/10.1175/1520-0469(1996)053:3468:ENMVEO>2.0.CO;2.

Bühler, Oliver (2014). *Waves and Mean Flows*. 2nd ed. Cambridge Monographs on Mechanics. Cambridge University Press. DOI: https://doi.org/10.1017/CBO9781107478701.

Burns, Keaton J., Geoffrey M. Vasil, Jeffrey S. Oishi, Daniel Lecoanet, and Benjamin P. Brown (Apr. 2020). "Dedalus: A flexible framework for numerical simulations with spectral methods". In: *Physical Review Research* 2.2, 023068, p. 023068. DOI: https://doi.org/10.1103/PhysRevResearch.2.023068.

Dunphy, Michael, Aurelien Ponte, Patrice Klein, and Sylvie Le Gentil (2017). "Low-mode internal tide propagation in a turbulent eddy field". In: *Journal of Physical Oceanography* 47.3. Ed. by Amer Meteorological Soc, pp. 649–665. DOI: https://doi.org/10.1175/JPO-D-16-0099.1.

Egbert, G. D. and S. Y. Erofeeva (2021). "An Approach to Empirical Mapping of Incoherent Internal Tides With Altimetry Data". In: *Geophysical Research Letters* 48.24, e2021GL095863. DOI: https://doi.org/10.1029/2021GL095863.

Godfrey, Michael D. (1965). "An exploratory study of the bi-spectrum of economic time series". In: *Journal of The Royal Statistical Society Series C-applied Statistics* 14, pp. 48–69.

Karban, U., E. Martini, A.V.G. Cavalieri, L. Lesshafft, and P. Jordan (2022). "Self-similar mechanisms in wall turbulence studied using resolvent analysis". In: *Journal of Fluid Mechanics* 939, A36. DOI: https://doi.org/10.1017/jfm.2022.225.

Loève (1955). *Probability theory*. Van Nostrand.

Lumley, J.L. (1967). "The Structure of Inhomogeneous Turbulent Flows". In: *Atmospheric Turbulence and Radio Wave Propagation*. Ed. by A.M. Yaglom and V.I. Tartarsky, pp. 166–177.

Ochoa, José, Julio Sheinenbaum, and Aleph Jiménez (2011). "Lateral Friction in Reduced-Gravity Models: Parameterizations Consistent with Energy Dissipation and Conservation of Angular Momentum". In: *Journal of Physical Oceanography* 41.10, pp. 1894–1901. DOI: https://doi.org/10.1175/2011JPO4599.1.

Ponte, A. L. and P. Klein (2015). "Incoherent signature of internal tides on sea level in idealized numerical simulations". In: *Geophysical Research Letters* 42.5, pp. 1520–1526. DOI: https://doi.org/10.1002/2014GL062583.

Richman, James G., Brian K. Arbic, Jay F. Shriver, E. Joseph Metzger, and Alan J. Wallcraft (2012). "Inferring dynamics from the wavenumber spectra of an eddying global ocean model with embedded tides". In: *Journal of Geophysical Research: Oceans* 117.C12. DOI: https://doi.org/10.1029/2012JC008364.

Savva, Miles and Jacques Vanneste (Dec. 2018). "Scattering of internal tides by barotropic quasigeostrophic flows". In: *Journal of Fluid Mechanics* 856, pp. 504–530. DOI: https://doi.org/10.1017/jfm.2018.694.

Schmidt, Oliver T. and Tim Colonius (2020). "Guide to Spectral Proper Orthogonal Decomposition". In: *AIAA Journal* 58.3, pp. 1023–1033. DOI: https://doi.org/10.2514/1.J058809.

Towne, Aaron, Oliver T. Schmidt, and Tim Colonius (2018). "Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis". In: *Journal of Fluid Mechanics* 847, pp. 821–867. DOI: https://doi.org/10.1017/jfm.2018.283.

Ward, Marshall L. and William K. Dewar (2010). "Scattering of gravity waves by potential vorticity in a shallow-water fluid". In: *Journal of Fluid Mechanics* 663, pp. 478–506. DOI: https://doi.org/10.1017/S0022112010003721.

Welch, P. (1967). "The use of fast Fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms". In: *IEEE Transactions on Audio and Electroacoustics* 15.2, pp. 70–73. DOI: https://doi.org/10.1109/TAU.1967.1161901.
**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_92.png" /></figure>
# Linear Wave Solutions of a Stochastic
Shallow Water Model

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

Check for
updates

Etienne Mémin, Long Li, Noé Lahaye, Gilles Tissot, and Bertrand Chapron

**Abstract** In this paper, we investigate the wave solutions of a stochastic rotating shallow water model. This approximate model provides an interesting simple description of the interplay between waves and random forcing ensuing either from the wind or coming as the feedback of the ocean on the atmosphere and leading in a very fast way to the selection of some wavelength. This interwoven, yet simple, mechanism explains the emergence of typical wavelength associated with near inertial waves. Ensemble-mean waves that are not in phase with the random forcing are damped at an exponential rate, whose magnitude depends on the random forcing variance. Geostrophic adjustment is also interpreted as a statistical homogenization process in which, in order to conserve potential vorticity, the small-scale component tends to align to the velocity fields to form a statistically homogeneous random field.

## 1 Introduction

Oceanic global circulation currents show a predominance of near inertial waves (NIW) in their spectrum. These waves result from the repeated forcing of atmospheric winds over the globe together with the influence of global earth rotation at Coriolis frequency. Besides, large-scale numerical representations of oceanic circulation require to introduce subgrid scale models that inescapably damp in the long run velocity fields and wave solutions. In such models high frequency waves are obviously completely smoothed out. But even large scale structures might be also affected if no particular attention is paid to the subgrid model design.

E. Mémin (✉) · L. Li · N. Lahaye · G. Tissot
ODYSSEY Team, Centre Inria de l'Université de Rennes, Rennes, France

IRMAR – UMR CNRS 6625, Rennes, France
e-mail: etienne.memin@inria.fr

B. Chapron
ODYSSEY Team, Centre Inria de l'Université de Rennes, Rennes, France
Laboratoire d'Océanographie Physique et Spatiale, Ifremer, Plouzané, France
In the last year, there has been an increasing effort to devise stochastic parameterizations for large-scale flows [2, 10, 11]. The motivations come from the failure of classical subgrid models to represent accurately the effect of the flow state variables at the unresolved scales and the necessity to provide reliable and computationally efficient models at the climatic scale. Uncertainty quantification, ensemble methods for forecasting and data assimilation are also prevailing, and the Bayes principle on which they are built leads naturally to considering stochastic dynamics.

Stochastic dynamics can hardly be devised on ad hoc grounds if one wants to provide generic and flexible systems. Control of the variance's growth and the respect of the physical properties of the underlying turbulent flows is of the utmost importance to define stable and physically relevant systems, in which the unresolved variables or the different physical and numerical approximations performed are faithfully represented. Two main schemes have been recently proposed to that purpose in the literature [12, 20]. The first one is a geometric framework relying on a Hamiltonian formulation, whereas the second one—referred to as modelling under location uncertainty (LU)—is based on Newton's principles and a stochastic formulation of the Reynolds transport theorem. Both schemes have been analyzed and numerically experimented for several geophysical models and configurations [1, 4, 5, 6, 17, 23, 25]. Surface waves and linear models have been proposed in the LU framework [8, 26, 27].

Nevertheless, the analysis of basic geostrophic mechanisms such as geostrophic adjustment or the form of basic wave solutions in linearized simple systems has not been investigated so far. This is the objective of this work. We will in particular analyze a stochastic rotating shallow water model recently proposed in [4, 16]. The focus will be on wave solutions of this system. We will show in particular that ensemble-mean waves ensuing from a random forcing at given frequencies are preserved whereas the others are very quickly damped. The stochastic model will allow us to propose a very simple mechanism leading to the emergence and conservation of NIW in atmosphere and ocean dynamics as a coupled self entrainment process.

The paper is organized as follows. In a first section we briefly recall the proposed stochastic model—rotating shallow water under location uncertainty (RSW-LU). We then identify a stationary solution of this stochastic system. We next analyze and illustrate the plane wave solutions associated to the linearized system under some specific noises. Finally, we provide a picture of the geostrophic adjustment process associated to the random system.

## 2 Review of RSW-LU

Let us first recall the stochastic transport operator introduced in the LU framework [1, 18, 17, 24]:
$$ \mathbb{D}_t \Theta := d_t \Theta + \left( \left( \mathbf{u} - \frac{1}{2} \nabla \cdot \mathbf{a} \right) dt + \sigma_t d\mathbf{B}_t \right) \cdot \nabla \Theta - \frac{1}{2} \nabla \cdot \left( \mathbf{a} \nabla \Theta \right) dt = 0, \quad (1a) $$

where the tracer, $\Theta$, is a stochastic process with an extensive property (e.g. temperature, salinity, buoyancy), $d_t \Theta(\mathbf{x}) := \Theta(\mathbf{x}, t + \delta t) - \Theta(\mathbf{x}, t)$ stands for the time-increment of $\Theta$ at a fixed point $\mathbf{x}$ with $\delta t$ an infinitesimal time variation, $\mathbf{u}$ denotes the time-smooth resolved velocity that is both spatially and temporally correlated, $\sigma_t d\mathbf{B}_t$ stands for the highly oscillating unresolved noise component, assumed in the present study to be divergence-free, with its quadratic variation [22] denoted by $\mathbf{a}$ and $\frac{1}{2} \nabla \cdot \mathbf{a}$ is the so-called Itô-Stokes drift [1] ensuing from the inhomogeneity of the noise. The mathematical definitions of the noise and its quadratic variation are given by

$$ \sigma_t d\mathbf{B}_t(\mathbf{x}) = \int_{\mathcal{D}} \check{\sigma}(\mathbf{x}, \mathbf{y}, t) d\mathbf{B}_t(\mathbf{y}) d\mathbf{y} = \sum_{n \in \mathbb{N}} \lambda_n^{1/2}(t) \xi_n(\mathbf{x}, t) d\beta_t^n, \quad (1b) $$

$$ a(\mathbf{x}, t) = \int_{\mathcal{D}} \check{\sigma}(\mathbf{x}, \mathbf{y}, t) \check{\sigma}^{\dagger}(\mathbf{y}, \mathbf{x}, t) d\mathbf{y} = \sum_{n \in \mathbb{N}} \lambda_n(t) (\xi_n \xi_n^{\dagger})(\mathbf{x}, t), \quad (1c) $$

where $\sigma_t$ is an integral operator defined on the Hilbert space $(L^2(\mathcal{D}))^d$ with a bounded spatial domain $\mathcal{D} \subset \mathbb{C}^d$ ($d = 2$ or 3), $\check{\sigma} = (\check{\sigma}_{ij})_{i,j=1,...,d}$ is a spatially and temporally bounded matrix kernel of $\sigma_t$, $\lambda_n$ and $\xi_n$ are the eigenfunctions and eigenvalues of the composite operator $\sigma_t \sigma_t^*$ ($\sigma_t^*$ denotes the adjoint of $\sigma_t$), $\bullet^\dagger$ stands for transpose-conjugate operation, $\mathbf{B}_t$ is a cylindrical Wiener process [22] and $\beta_t^n$ are independent (one-dimensional) standard Brownian motions.

Under the stochastic transport notations (1), the governing equations of the energy-preserved RSW-LU system [4] read

(Conservation of momentum)

$$ \mathbb{D}_t \mathbf{u} + f_0 \mathbf{u}^{\perp} dt = -g \nabla \eta dt, \quad (2a) $$

(Conservation of mass)

$$ \mathbb{D}_t h + h \nabla \cdot \mathbf{u} dt = 0, \quad (2b) $$

(Incompressible constraints)

$$ \sigma_t d\mathbf{B}_t = \nabla^{\perp} \varphi d\mathbf{B}_t, \quad \nabla \cdot \nabla \cdot \mathbf{a} = 0, \quad (2c) $$

(Conservation of energy)

$$ d_t \int_{\mathcal{D}} \frac{1}{2} \rho (h |\mathbf{u}|^2 + g h^2) d\mathbf{x} = 0, \quad (2d) $$

where $\mathbf{u} = [u, v]^T$ denotes the two-dimensional horizontal velocity with $\mathbf{u}^{\perp} = [-v, u]^T$, $\nabla = [\partial_x, \partial_y]^T$ stands for the horizontal gradient with $\nabla^{\perp} = [-\partial_y, \partial_x]^T$,
$h(\mathbf{x}, t) = H + \eta(\mathbf{x}, t)$ is the water thickness with $H$ the flat bottom height and $\eta$ the free-surface position, $f_0$ is a constant approximation of the Coriolis parameter, $g$ is the gravity constant, $\rho$ is the water density and $\varphi dB_t$ denotes a random (scalar) stream function defined in a similar way as in (1b). As shown in [4], the incompressible conditions (2c) for both noise and Itô-Stokes drift ensure the path-wise conservation of the total energy (2d). We remark that the analytical properties for the RSW-LU system have been investigated in [16] and some numerical applications of (2) have been performed in [4].

### 3 Stationary Solution

We focus now on stationary solutions associated to the previous system. To that end, neglecting the time increments $d_t u$ and $d_t h$ in (2a) and (2b) and recalling that $\nabla \eta = \nabla h$ due to the flat bottom assumption, we obtain

$$ \left( \left( u - \frac{1}{2} \nabla \cdot a \right) \cdot \nabla u - \frac{1}{2} \nabla \cdot \left( a \nabla u \right) + f_0 u^{\perp} + g \nabla h \right) dt + \sigma_t dB_t \cdot \nabla u = 0, \quad (3a) $$

$$ \left( \left( u - \frac{1}{2} \nabla \cdot a \right) \cdot \nabla h - \frac{1}{2} \nabla \cdot \left( a \nabla h \right) + h \nabla \cdot u \right) dt + \sigma_t dB_t \cdot \nabla h = 0, \quad (3b) $$

$$ \sigma_t dB_t = \nabla^{\perp} \varphi dB_t, \quad a = \nabla^{\perp} \varphi (\nabla^{\perp} \varphi)^T. \quad (3c) $$

From Bichteler-Dellacherie decomposition of a semi-martingale, the martingale terms (i.e. the Brownian terms) and the finite variation terms (i.e. the differentiable terms) can be safely separated. Decomposing in such way the mass equation (3b) and considering the corresponding martingale part leads to $\sigma_t dB_t \cdot \nabla h = 0$. The random surface height gradient is aligned with the large-scale surface height gradient

$$ \sigma_t dB_t = \alpha \nabla^{\perp} h, \quad (4) $$

The noise being divergence free yields that $\alpha$ is constant along the level sets of $h$:

$$ \nabla \alpha \cdot \nabla^{\perp} h = 0. \quad (5) $$

It is hence necessary of the form $\alpha = \phi(h) + C$, with $C$ a constant, and $\phi$ a differentiable function. For such a noise the variance tensor is given by

$$ a = \alpha^2 \nabla^{\perp} h \nabla^{\perp} h^T, \quad (6) $$

while the ISD reads
$$ \nabla \cdot a = \alpha^2 \nabla \cdot (\nabla^{\perp} h \nabla^{\perp} h^T) = \alpha^2 (\nabla^{\perp} h \cdot \nabla) \nabla^{\perp} h, \quad (7) $$

and we notice that $a\nabla h = 0$. The mass equation then boils down to

$$ (u - \frac{1}{2} \nabla \cdot a) \cdot \nabla h + h \nabla \cdot u = 0. \quad (8a) $$

In the same way, we get from a semi-martingale decomposition of (3a) that

$$ \sigma_t dB_t \cdot \nabla u = 0, \quad (8b) $$

which yields

$$ (u - \frac{1}{2} \nabla \cdot a) \cdot \nabla u + f_0 u^{\perp} + g \nabla h = 0. \quad (8c) $$

At this point it is interesting to interpret the random streamfunction in terms of nondimensional units to infer the importance of the Itô-Stokes drift (ISD) in comparison to geostrophic flows. To that end, we assume that $\varphi dB_t \sim \sqrt{\epsilon T} (g H/f_0) \sqrt{T}$, where $T$ denotes the characteristic time scale and $\epsilon T$ stands for a small-scale decorrelation time with $\epsilon$ the strength of uncertainty (the greater $\epsilon$ the stronger the noise is). From definition (3c), the noise's quadratic variation scales then as $a \sim \epsilon T (g H)^2 / (f_0 L)^2$ with $L$ the characteristic length scale. As a consequence, the ratio of the ISD advection term with the geostrophic gradient pressure scales as

$$ \frac{(\nabla \cdot a) \cdot \nabla u}{g \nabla h} \sim \epsilon \frac{g H}{f_0^2 L^2} = \epsilon \frac{L_d^2}{L^2} = \epsilon B_u, \quad (9) $$

where $L_d = \sqrt{gH}/f_0$ is the Rossby deformation radius, and $B_u$ denotes the Burger number, which stands for the ratio between vertical density stratification and the earth's rotation in the horizontal ($B_u = (NH/\Omega L)^2 = (L_d/L)^2$). If the ISD advection has the same importance as the gradient pressure term, then we must have $\sqrt{\epsilon} = L/L_d$. This means that when the scale of motions is greater than the deformation radius, the rotation will dominate and the noise must be important to have the ISD playing a role. At the opposite, when the scale of motions is smaller than the Rossby radius, the small-scale flow component does not need to be important to be as significant as the gradient pressure term. From this point of view, the deformation radius can be interpreted as the limiting scale under which the statistically modified advection due to the inhomogeneity of the small-scale component plays a role.

In the following, we will assume to be at a scale much larger than the deformation radius, so that the action of the ISD becomes negligible. In that case, the stationary system finally simplifies as

$$ (u \cdot \nabla)u + f_0 u^{\perp} + g \nabla h = \nabla \cdot (uh) = 0. \quad (10) $$
This system is structurally the same as the nonlinear stationary system, at the exception of a scaling constraint on the ISD, which makes negligible the nonlinear advection term. As a matter of fact, noticing that for $\mathbf{u} \propto \nabla^{\perp}h$, the advection term corresponds to the ISD (7), it follows that the geostrophic balance, $\mathbf{u} = -(g/f_0)\nabla^{\perp}h$, is a stationary solution of such system for a null ISD. We show below it is the only non-trivial stationary solution of such system.

Let us write the velocity as a superposition of the geostrophic component and an ageostrophic component:

$$ \mathbf{u} = -\frac{g}{f_0} \nabla^{\perp}h + \mathbf{v}, \qquad (11) $$

where $\mathbf{v}$ is defined through the Helmholtz decomposition from a potential function, $\Phi$, and a stream function $\Psi$ that both depend on the surface elevation:

$$ \mathbf{v} = \nabla^{\perp}\Psi(h) + \nabla\Phi(h), \qquad (12) $$

$$ = \nabla^{\perp}h \Psi'(h) + \nabla h \Phi'(h). \qquad (13) $$

As $\nabla^{\perp}h$ belongs to the null space of the velocity gradient (8b), $\nabla h$ either belongs also to the null space or it is an eigenvector of the velocity gradient tensor ($\nabla \mathbf{u}$). From the momentum equation (10) we have

$$ (\nabla h \Phi'(h)) \cdot \nabla \mathbf{v} = -f \mathbf{v}^{\perp}. $$

Then if $\nabla h$ belongs to the null space of the velocity gradient, $\mathbf{v}$ directly cancels out. If is an eigenvector of the velocity gradient with eigenvalue $\lambda$, the above equation reads

$$ \lambda \nabla h \Phi'(h) = -f(\nabla h \Psi'(h) + \nabla^{\perp}h \Phi'(h)), $$

which implies $\Phi'(h) = \Psi'(h) = 0$ and hence $\mathbf{v} = 0$.

Physically, we see hence that the considered nonlinear system admits a stationary solution for a negligible ISD. Let us note that no such constraint is available in the deterministic setting. In the following we will interpret this stationary solution in the linearized stochastic shallow water system in terms of wave solutions and geostrophic adjustment.

# 4 Stochastic Rotating Shallow Water Waves

In order to look at the different wave solutions associated to the stochastic shallow water system (2), we proceed in the same way as in the deterministic case [19, 28]. In particular, we assume that the noise's structure in (1b) is independent of the resolved
prognostic variable $\mathbf{u}$. In that case, the associated linearized system of (2) reads

$$d_t \mathbf{u} + (f_0 \mathbf{u}^{\perp} + g \nabla \eta - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \mathbf{u})) dt + \sigma_t d\mathbf{B}_t \cdot \nabla \mathbf{u} = 0, \quad (14a)$$

$$d_t \eta + \left( H \nabla \cdot \mathbf{u} - \frac{1}{2} \nabla \cdot (\mathbf{a} \nabla \eta) \right) dt + \sigma_t d\mathbf{B}_t \cdot \nabla \eta = 0, \quad (14b)$$

$$\nabla \cdot \sigma_t d\mathbf{B}_t = \nabla \cdot \mathbf{a} = 0, \quad (14c)$$

$$d_t \int_D \frac{1}{2} \rho (H |\mathbf{u}|^2 + g (H + \eta)^2) d\mathbf{x} = 0. \quad (14d)$$

It can be checked that this system conserves the total energy (14d) in the same way as the initial nonlinear stochastic system (2) does. For noise defined through Hilbert-Schmidt correlation tensor this system admits a mild solution [22]. Existence of strong solution could also be shown from the nonlinear system [16]. In order to build some simple analytical solutions of this linearized system and to better understand the physical behaviours of the resulting waves, only very specific noise models will be considered in the following. We first build the ensemble-mean wave solutions under homogeneous noise, then investigate the path-wise solutions under constant noise and under homogeneous noise but with very smooth structures.

## 4.1 Ensemble-Mean Waves Under Homogeneous Noise

Let us first recall the definition of the homogeneous and incompressible noise:

$$\begin{aligned} \sigma_t(\mathbf{x}) d\mathbf{B}_t &= \int_D \nabla^\perp \varphi(\mathbf{x} - \mathbf{y}) d\mathbf{B}_t(\mathbf{y}) d\mathbf{y} \\ &= \sum_m i \mathbf{k}_m^\perp \hat{\varphi}(\mathbf{k}_m) d\beta_t(\mathbf{k}_m) \exp(ik_m \cdot \mathbf{x}), \end{aligned} \quad (15a)$$

$$\mathbf{a} = \sum_m |\hat{\varphi}(\mathbf{k}_m)|^2 \mathbf{k}_m^\perp (\mathbf{k}_m^\perp)^T, \quad (15b)$$

where $i$ denotes the imaginary unit, $\mathbf{k}_m = [k_m, \ell_m]^T \in \mathbb{R}^2$ is the $m$-th wavenumber vector, $\hat{\varphi}$ stands for Fourier transform (in space) coefficient and $\beta_t \in \mathbb{C}$ are independent Brownian motions satisfying $\beta_t(-\mathbf{k}_m) = \bar{\beta}_t(\mathbf{k}_m)$ with Re{$\beta_t$} and Im{$\beta_t$} be independent. This noise is homogeneous, and thus associated to a constant matrix $\mathbf{a}$. Its ISD is null and fits naturally the condition on the stationary solution found in the previous section.

For sake of simplicity, we assume hereafter that the noise is defined by only one Fourier mode (associated to a wavenumber $\mathbf{k}_\sigma$) combined with its complex