# Physically constrained covariance inflation from location uncertainty

Yicun Zhen, Valentin Resseguier, Bertrand Chapron

► To cite this version:

Yicun Zhen, Valentin Resseguier, Bertrand Chapron. Physically constrained covariance inflation from location uncertainty. Nonlinear Processes in Geophysics, 2023, 30 (2), pp.237-251. 10.5194/npg-30-237-2023. hal-04204126

HAL Id: hal-04204126

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

Submitted on 13 Sep 2023

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

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

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

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

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

# Physically constrained covariance inflation from location uncertainty

Yicun Zhen<fnref n="1" />, Valentin Resseguier<fnref n="2" />,<fnref n="3" />, and Bertrand Chapron<fnref n="4" />

**Correspondence:** Yicun Zhen (zhenyicun@proton.me)

Received: 6 March 2023 – Discussion started: 9 March 2023

Revised: 13 May 2023 – Accepted: 23 May 2023 – Published: 29 June 2023

**Abstract.** Motivated by the concept of “location uncertainty”, initially introduced in Mémin (2014), a scheme is sought to perturb the “location” of a state variable at every forecast time step. Further considering Brenier’s theorem (Brenier, 1991), asserting that the difference of two positive density fields on the same domain can be represented by a transportation map, we demonstrate that the perturbations consistently define a stochastic partial differential equation (SPDE) from the original PDE. It ensues that certain quantities, up to the user, are conserved at every time step. Remarkably, derivations following both the SALT (stochastic advection by Lie transport; Holm, 2015) and LU (location uncertainty; Mémin, 2014; Resseguier et al., 2017a) settings can be recovered from this perturbation scheme. Still, it offers broader applicability since it does not explicitly rely on Lagrangian mechanics or Newton’s laws of force. For illustration, a stochastic version of the thermal shallow water equation is presented.

## 1 Introduction

Data assimilation is meant to extract information from measurements to improve the state estimate. Kalman-filter-based and particle-filter-based methods are now commonly used for academic studies and operational forecasts. For both methods, the estimate of a state variable and the uncertainty quantification of the estimate of a state variable are repeated at each data assimilation cycle. In the classical Kalman filter, this uncertainty is represented by a covariance matrix. In Monte Carlo-based methods (i.e., the ensemble Kalman filters and particle filters, etc.), it is represented by the spread of the ensemble members or particles. The uncertainty of the state estimate is further part of the input for the next data assimilation cycle. Frequently observed, the uncertainty can be underestimated in nonlinear numerical experiments when there is no model noise (Schlee et al., 1966; Harlim and Majačka, 2010; Franzke et al., 2015). As a consequence, the state estimate in the subsequent time steps may not be efficiently adjusted by the physical measurements: the system is over-

confident about its current state estimate. This phenomenon is usually referred to as filter divergence, possibly associated with the “curse of dimensionality” (see for instance Daum and Huang, 2003).

To address the latter issue, “covariance localization” has been developed for both Kalman-filter-based methods and particle filters (Houtekamer and Mitchell, 2001; Poterjoy, 2016). To further mitigate filter divergence, a practical strategy is to inflate the uncertainty estimate at each forecast time step or each data assimilation cycle (Anderson, 2007; Tibshirani and Knight, 1999; Li et al., 2009; Kotsuki et al., 2017; Ying and Zhang, 2015; Miyoshi, 2011; Raanes et al., 2019; Zhen and Harlim, 2015). For geophysical applications, the uncertainty is then often inflated by rescaling the ensemble covariance in order to match bias and variance. A natural alternative is the addition of noise in the dynamical equations.

In the context of ensemble-/particle-based methods, the uncertainty is usually inflated by artificially perturbing each ensemble member/particle. We refer the reader to

<footnote n="1">¹College of Oceanography, Hohai University, Nanjing, China</footnote>

<footnote n="2">²LAB SCALIAN DS, Rennes, France</footnote>

<footnote n="3">³INRAE, OPAALE, Rennes, France</footnote>

<footnote n="4">⁴Laboratoire d'Océanographie Physique et Spatiale, Ifremer, Plouzaé, France</footnote>
Resseguier et al. (2021) for a review on the subject. It is then a natural question to ask whether there is a mathematical principle to guide this uncertainty inflation.

In the fluid dynamics community, random forcings are not introduced for inflation but to mimic the intermittent backscattering of energy from small scales toward large scales. Among those approaches, we may mention the stochastic Lagrangian models (Pope, 1994) and the Eulerian Gaussian backscatterings of the EDQNM (eddy damped quasi-normal Markovian; Orszag, 1970; Leith, 1971) model. Additive noise models, like the linear inverse models (Peland and Sardeshmukh, 1995), have then also been proposed for filtering purposes and have been thoroughly reviewed by Tandeo et al. (2020). Most methods mainly focus on comparing the estimated uncertainty and the statistics of the innovation process but ignore other mathematical/physical aspects (for instance, conservation laws). Other empirical approaches, referred to as the SPPT (stochastically perturbed parametrization tendency; Buizza et al., 1999) scheme and the SKEBS (stochastic kinetic energy backscatter; Berner et al., 2009) scheme, introduce multiplicative noise, with success in operational weather and climate forecast centers (Franzke et al., 2015). Still many drawbacks have been reported, above all violations of conservation laws (Reynolds et al., 2016; Leutbechner et al., 2016). Recently, the operational ocean circulation model NEMO has also been randomized (e.g. Leroux et al., 2022) but, again, without consideration of conservation.

Several authors proposed schemes specifically to enforce energy conservation or at least a given energy budget (e.g. Sapsis and Majda, 2013; Gugole and Franzke, 2019; Resseguier et al., 2021). To better constrain non-Gaussian schemes, many authors rely on physics and possibly on time-scale separation. Introduced by Hasselmann (1976), it is generally associated with the rigorous theories of averaging and homogenization. Majda et al. (1999) decomposed the state variable into slow-varying modes $x_j$ and fast-varying modes $y_j$. The authors demonstrated that the interaction term between $x_j$ and $y_j$, in the equation for $x_j$, can be modeled as a stochastic process solely in terms of $x_j$'s, as the ratio of the timescales of $x_j$ and $y_j$ tends to 0. Nevertheless, homogenization methods, like Majda et al. (1999), may also lead to violation of energy conservation, even though some workarounds exist (Frank and Gottwald, 2013; Jain et al., 2014).

In Brzeźniak et al. (1991), later modified in Mikulevicius and Rozovskii (2004), Flandoli (2011), Mémin (2014), and Resseguier et al. (2017a, 2021), preservation of kinetic energy is specifically emphasized. The true velocity of an incompressible flow is decomposed into a regular component and a turbulent one and the latter modeled by a stochastic noise. Mikulevicius and Rozovskii (2004) and Mémin (2014) further derived stochastic Navier–Stokes equations. For these two approaches, the large-scale advecting velocity differs, induced by different regularization of Newton' second law.

Following another path, considering Hamilton’s principle with a stochastic advection constraint on Lagrangian fluid trajectories, Holm (2015) also proposed a consistent stochastic setting, i.e., stochastic advection by Lie transport (SALT). In particular, this derivation preserves Kelvin’s circulation. Similarities and differences between these different stochastic frameworks are discussed in Resseguier et al. (2020).

From another perspective, the classical optimal transport theory suggests that the difference of two smooth positive density fields ($ρ_1$ and $ρ_2$) on a bounded domain $Ω$ can be described by a transportation map: $T: Ω → Ω$. More specifically, there exists a diffeomorphism $T$ of $Ω$ to transform $ρ_1$ to $ρ_2$ under the diffeomorphism $T$ with a minimal cost. Broadly speaking, $T$ can be interpreted as how much $ρ_2$ differs from $ρ_1$, and $T$ operates as a location correction. Indeed, starting from the same initial condition $ρ(t)$, suppose that $ρ_1 = ρ^{\text{model}}(t + Δt)$ is the model forecast and $ρ_2 = ρ(t + Δt)$ is the true forecast. The additional uncertainty of $ρ_1$ due to model error can then be represented by a random $T$. It further suggests that the inflation of uncertainty can be achieved by casting a random $T$ on each ensemble member/particle.

Motivated by such an optimal transport perspective and the concept of “location uncertainty”, proposed in Mémin (2014), a new strategy can thus seek to design a well-constrained “location perturbation” of the state variable.

Specifically, the idea of covariance inflation can be informally generalized to physical fields that are not always positive, i.e., physical fields other than the density field. Mathematically, a density field $ρ$ is naturally associated with a differential $n$ form $θ_ρ$, where $n = \text{dim}Ω$. The statement “$ρ_1$ transforms to $ρ_2$ under the diffeomorphism $T$” is equivalent to the mathematical relation $θ_{ρ_1} = T^*θ_{ρ_2}$, where $T^*$, acting on all differential forms, is the pull-back operator induced by $T$ or, equivalently, $θ_{ρ_2} = (T^{-1})^*θ_{ρ_1}$. Therefore, a random $T$ (or equivalently, $T^{-1}$) could induce a perturbation of any differential $k$ form.

To implement a physically constrained perturbation scheme, the state variable $S$ under consideration must then be associated with some differential form $θ$, i.e., construct a one-to-one correspondence between snapshots of $S$ and snapshots of $θ$. Note that this can be generalized to other types of tensor fields.

It will be demonstrated (Sect. 5) that it is indeed sometimes helpful to choose $θ$ to be a contravariant tensor field other than differential forms. Yet, it must be stressed that associating the state variable $S$ with a differential form $θ$ is a key important step.

Correspondingly, at each forecast time step, the covariance inflation should follow four steps.

- Step 1: find $θ(t)$ based on $S(t)$.

- Step 2: construct a random diffeomorphism $T: Ω → Ω$.

- Step 3: replace $θ(t)$ with $T^*θ(t)$ and calculate $S(t)$ based on the new value of $θ(t)$.
- Step 4: calculate the forecast $S(t+\Delta t)$ based on the new value of $S(t)$.

Associating $S$ with different $\theta$ shall then be constrained by different conservation laws for the perturbation scheme. More precisely, certain physical quantities are conserved in step 3, no matter how $T$ is constructed or realized in step 2. We emphasize that the conservation law of the perturbation scheme merely depends on the choice of $\theta$ but is independent of the dynamics of the original deterministic system. A resulting stochastic partial differential equation (SPDE) will conserve a given quantity only if both the perturbation scheme and the original deterministic system conserve that quantity. We also remark that this scheme can not conserve all the physical quantities at the same time unless additional constraints upon the parameters are imposed. Hence the users must choose by themselves which physical quantity to conserve.

In summary, the key perspective of this paper is that the displacement vector field of physical state variables should be determined by the tensor fields associated with the physical fields. The advantage of this perspective is that certain physical quantities can be conserved while applying a displacement vector field to transfer the original physical field. A direct application of this perspective is the physically constrained covariance inflation scheme proposed in this paper. When the tensor fields are positive $n$ forms on a bounded domain that have the same total mass, Brenier's theorem shows that the "optimal" displacement vector field exists and is unique, for a given cost function. In this case, the optimality of displacement vector field is well-defined. In other cases, the issue of "optimality" together with the existence and uniqueness of the optimal displacement vector field need to be carefully explored. We reserve this for future study.

This paper is organized as follows. Section 2 is a brief introduction of optimal transport theory. In Sect. 3 we present the perturbation scheme in detail, including the motivation, the specific techniques in derivation, and several examples. In Sect. 4, the resulting perturbation scheme is then compared with the stochastic advection by Lie transport (SALT) equations (Holm, 2015) and the location uncertainty (LU) equations (Mémin, 2014). For properly chosen $\theta$ and $T_t$, it is demonstrated that both SALT and LU settings are recovered within the proposed framework. To illustrate our purpose, a stochastic version of the thermal shallow water equation is then derived in Sect. 5. A final conclusion and discussion is given in Sect. 6.

The conventions of notation are as follows:

- The letter $i$ only refer to the $i$th independent Brownian motion. The letters $p, q, j$, and $k$ refer to the components if $p, q, j$, and $k$ are upper indices.

- For Einstein's convention on summation (applies to all indices except $i$ and $j$), if $p$ is shown in both upper and lower indices, then the summation over $p$ automatically applies.

- Summation over $i, j$, and $p$ automatically applies in all equations. For instance, $e_i$ refers to $\sum_i e_i$, and $y_j$ refers to $\sum_j y_j$.

## 2 Monge's formulation of optimal transport problem and Brenier's answer

Hereafter we briefly summarize some necessary concepts and results in optimal transport theory. Let $\Omega$ be a bounded domain in a $n$-dimensional Euclidean space.

**Definition 2.0.1** (Monge's optimal transport problem). Given the cost function $c(x, y) \ge 0$ and probability measures $\mu, \nu \in \mathcal{P}(\Omega)$,

$$ \text{minimize} \quad \mathbb{M}(T) = \int_{\Omega} c(x, T(x)) d\mu(x) \quad (1) $$

over $\mu$ measurable maps $T: \Omega \to \Omega$ subject to $\nu = T_{\text{no.}\mu}$.

Here the probability measures $\mu$ and $\nu$ are interpreted as mass distributions with total mass equal to 1. The map $T$ is called a transport plan which moves the mass $d\mu(x)$ at location $x$ to location $T(x)$, with the cost $c(x, T(x))$ per unit of mass. Therefore the quantity $\mathbb{M}(T)$ is the total cost of the transport plan $T$. The constraint $\nu = T_{\text{no.}\mu}$ is interpreted as $T$ transporting the mass distribution $\mu$ to the mass distribution $\nu$. In the case that $T$ is a diffeomorphism and that both $\nu$ and $\mu$ have smooth densities, i.e., assuming that $dv(x) = f(x)d^n x$ and $d\mu(x) = g(x)d^n x$ for some smooth functions $f, g$ on $\Omega$,

$$ \nu = T_{\text{no.}\mu} \iff g(x) = f(T(x))|J_T(x)|, \quad (2) $$

where $J_T(x)$ refers to the Jacobian matrix of $T$ at $x$. If we associate $\nu$ and $\mu$ with differential $n$ forms $\theta_\nu = f dx^1 \wedge \cdots \wedge dx^n$ and $\theta_\mu = g dx^1 \wedge \cdots \wedge dx^n$, then

$$ \nu = T_{\text{no.}\mu} \iff \theta_\mu = T^* \theta_\nu. \quad (3) $$

Brenier (1991) proved the existence and uniqueness of the solution to Monge's optimal transport problem for $c(x, y) = |x - y|^2$. To better illustrate how optimal transport theory motivates us, we consider the following simplified version of Brenier's theorem.

**Theorem 1** (Brenier, simplified version). Let $\mu$ and $\nu$ be measures with bounded smooth density on a bounded domain $\Omega \subset \mathbb{R}^n$. Let $c(x, y) = |x - y|^2$. Then there is a convex function $\phi: \Omega \to \mathbb{R}$, such that $(\nabla\phi)_{\text{no.}\mu} = \nu$. And $\nabla\phi: x \to x + \nabla\phi|_x$, defined $\mu$-almost everywhere, is the unique solution to Monge's optimal transport problem.

The convexity of $\phi$ implies that the map $\nabla\phi$ is one to one. Broadly speaking, Brenier's theorem implies that the difference of two density fields can be represented by a transportation map $T$.
### 3 The perturbation scheme

Consider a compressible flow on a bounded domain $\Omega$. Let $\rho$ denote the density field. Let $\rho^{\text{model}}(t + \Delta t)$ and $\rho^{\text{true}}(t + \Delta t)$ be the model forecast and the true forecast starting from the same density field at time $t$. If we assume that the model forecast and the truth have the same total mass, Brenier's theorem says that there exists a diffeomorphism $T: \Omega \to \Omega$, so that

$$ \rho^{\text{true}}(x, t + \Delta t) = \rho^{\text{model}}(T(x), t + \Delta t)J_T(x). \quad (4) $$

Note that the transportation $T$ hereinafter is equivalent to the mapping $T^{-1}$ used in the Introduction. Equation (4) can further be written in terms of differential form. Let $\theta_{\rho} = \rho dx^1 \wedge \dots \wedge dx^n$, then Eq. (4) is equivalent to

$$ T^* \theta_{\rho}^{\text{model}}(t + \Delta t) = \theta_{\rho}^{\text{true}}(t + \Delta t). \quad (5) $$

For general differential forms $\theta$, it is unclear whether a diffeomorphism $T$ always exists that satisfies Eq. (5). However, Eq. (5) provides us with a tool for covariance inflation by constructing a random $T$ at every infinitesimal time step.

At each time step we construct a small perturbation $T$:

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

where $a(t, x)$, $e_i(t, x) \in \mathbb{R}^n$, $\Delta \eta_i(t) \sim \mathcal{N}(0, \Delta t)$ is a random number. Essentially, $T_t(x) - x$ can be interpreted as a “location error” caused by the model error. In Eq. (6), $a(t, x)\Delta t$ refers to a systematic location error, and $e_i \Delta \eta_i$ refers to a random location error.

Stated in the Introduction, the state variable $S$ must first be associated with a differential form $\theta$. Then at every time step, $T_t$ induces a perturbation of $\theta(t)$ by $\theta(t) \to T_t^*\theta(t)$. It hence induces a perturbation of the state variable $S(t)$. A forecast is then performed based on the perturbed state. Consequently, this perturbation scheme derives an SPDE from the original PDE.

This procedure can also be generalized to other types of tensor fields. We refer to Chern et al. (1999) for a rigorous definition of the tensor fields and the wedge algebra. For instance, we may choose $\theta = \rho \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}$, where $\{\frac{\partial}{\partial x^i}\}_{i \le n}$ forms a global basis of the tangent field. Then $T_t$ induces a perturbation of $\theta$ by $\theta(t) \to T_{t*}\theta$, where $T_{t*}$ is the push-forward operator induced by $T_t$. In Sect. 5, such a generalization is found useful in the example of thermal shallow water equation.

**Remark 1.** When $\theta$ is a mixture of covariant and contravariant tensor fields, the perturbation scheme is slightly more complicated. Assume that $T_t: \Omega_1 \to \Omega_2$ is a diffeomorphism, and $\theta = v \otimes \omega$, where $v$ and $\omega$ are contravariant or covariant tensor fields respectively on $\Omega_2$. Then $T_t^*\omega$ is a covariant tensor field on $\Omega_1$. However, $T_t$ can not directly induce a contravariant tensor field on $\Omega_1$. In order to get a tensor field on $\Omega_1$, we consider $T_t^{-1}: \Omega_2 \to \Omega_1$ and apply the push-forward operator on $v$. In sum, we may define the

perturbation to be

$$ \theta(t) \to ((T_t^{-1})_* v) \otimes (T_t^* \omega). \quad (7) $$

Appendix A derives the expression of $T_t^{-1}$ directly from the expression of $T_t$.

### 3.1 Calculation of $T_t^*\theta$ (or $T_{t*}\theta$)

A rigorous mathematical definition and calculation of $T_t$ and $T_t^*$ should be given in terms of stochastic flows of diffeomorphisms and its Lie derivatives. A brief discussion of the relationship between $T_t^*$ and the Lie derivative is given in Sect. 4.1. We further refer to Leon (2021) for a detailed definition of the Lie derivative. Yet, to rapidly assess $T_t^*\theta$ (or $T_{t*}\theta$), a Taylor expansion and Itô's lemma can be used.

Given coordinates $(x^1, \dots, x^n)$, when $\theta$ is a differential $k$ form, it can be written as

$$ \theta = \sum_{i_1 < \dots < i_k} f^{i_1, \dots, i_k} dx^{i_1} \wedge \dots \wedge dx^{i_k}. \quad (8) $$

Then

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

Given in Appendix B, a Taylor expansion and Itô's lemma are applied to expand $T_t^*\theta$, leading us to compactly write

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

for some differential $k$ forms $\mathcal{M}(\theta)$ and $\mathcal{N}_i(\theta)$. Hereafter, several examples of $T_t^*\theta$ are presented.

The full derivation of these examples is skipped. We further express all the terms in coordinates. For instance, we replace $\langle \nabla f \rangle a$ with $a^j \partial_{x^j} f$, where, by convention of notation, $a^j \partial_{x^j} f = \sum_j a^j \frac{\partial f}{\partial x^j}$. Similarly, $e_i^\top H_f e_i$ is replaced with $e_i^p e_i^q \partial_{x^p} \partial_{x^q} f$.

**Remark 2.** When $\theta = f \frac{\partial}{\partial x^{i_1}} \wedge \dots \wedge \frac{\partial}{\partial x^{i_k}}$ is a contravariant tensor field,

$$ T_{t*} \theta = f(T_t^{-1}(x)) T_{t*} \left( \frac{\partial}{\partial x^{i_1}} \wedge \dots \wedge \frac{\partial}{\partial x^{i_k}} \right). \quad (11) $$

The formula for $T_t^{-1}$ is derived in Appendix A. Then the expression of $f(T_t^{-1}(x))$, $T_{t*} \frac{\partial}{\partial x^{i_1}} \wedge \dots \wedge \frac{\partial}{\partial x^{i_k}}$ and $T_{t*}\theta$ can be derived step by step in a similar way to that in the Appendix B.

**Example 3.1.1** When $\theta = f$ is a function (differential 0 form),

$$ (T_t^* \theta) = f + (a^j \partial_{x^j} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f) \Delta t \\ + e_i^p \partial_{x^p} f \Delta \eta_i. \quad (12) $$

**Example 3.1.2** When $\theta = dx^1 \wedge dx^2 \wedge \dots \wedge dx^n$,

$$ T_t^* \theta = \{1 + (\partial_{x^p} a^p + \frac{1}{2} J_i) \Delta t + \partial_{x^p} e_i^p \Delta \eta_i\} \theta, \quad (13) $$
where $J_i = \partial_{x^p} e_i^p \partial_{x^q} e_i^q - \partial_{x^p} e_i^q \partial_{x^q} e_i^p$.

**Example 3.1.3** When $\theta = f \, dx^1 \wedge \dots \wedge dx^n$,

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

**Example 3.1.4** When $\theta = f^j \, dx^j$ (note that by the convention of notation, $f^j \, dx^j = \sum_{j=1}^n f^j \, dx^j$),

$$ T_t^* \theta = \{ f^j + (a^p \partial_{x^p} f^j + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f^j + \partial_{x^j} a^p f^p \\ + \partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p) \Delta t \\ + (e_i^p \partial_{x^p} f^j + \partial_{x^j} e_i^p f^p) \Delta \eta_i \} dx^j. \quad (15) $$

**Example 3.1.5** When $\theta = f \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}$,

$$ T_{t*} \theta = \{ f + ((\partial_{x^p} a^p + \frac{1}{2} J_i) f + (- (a^p + e_i^p \partial_{x^q} e_i^q) \\ + \partial_{x^q} e_i^p e_i^q) \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f) \Delta t \\ + (\partial_{x^p} e_i^p f - e_i^p \partial_{x^p} f) \Delta \eta_i \} \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}. \quad (16) $$

## 3.2 Derivation of the SPDE

Suppose S is the full state variable of the dynamical system:

$$ \frac{\partial S}{\partial t} = g(S). \quad (17) $$

Let f be a component or a collection of components of S. We then associate f with a differential form $\theta$ in the perturbation scheme; i.e., there is an invertible map $\mathcal{F}$ that maps the space of f to the space of $\theta$, such that $\mathcal{F}(f) = \theta$. In the examples in this paper, the corresponding map $\mathcal{F}$ is obvious. When f refers to a scalar quantity on the domain. We can choose to associate a differential n form $\theta = f \, dx^1 \wedge \dots \wedge dx^n$ with f, as in Example 3.1.3, an n vector $\theta = f \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}$ as in Example 3.1.5, or a differential 0 form (function) $\theta = f$. When f refers to a vector-valued function $f = (f^1, \dots, f^n)$, we can associate the differential 1 form $f = f^1 \, dx^1 + \dots + f^n \, dx^n$ with f. It is not hard to see that $\mathcal{F}$ is obvious once the type of tensor field is chosen. Suppose the propagation equation for f is

$$ d f = g^f(S) dt. \quad (18) $$

This implies a propagation equation for $\theta$:

$$ d\theta = g^\theta(S)dt. \quad (19) $$

The discrete-time perturbed forecast at each time step consists of the following two steps:

$$ \tilde{\theta}(t + \Delta t) = \theta(t) + g^\theta(S(t)) \Delta t \quad (20) $$

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

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

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

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

Therefore,

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

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

$$ d\theta = g^\theta(S)dt + \mathcal{M}(\theta)dt + \mathcal{N}_i(\theta)\Delta\eta_i. \quad (24) $$

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

$$ d f = g^f(S) dt + \mathcal{M}^f(f) dt + \mathcal{N}_i^f(f) d\eta_i. \quad (25) $$

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

$$ d_s f := \mathcal{M}^f(f) dt + \mathcal{N}_i^f(f) d\eta_i. \quad (26) $$

Then Eq. (25) can be written as

$$ d f = g^f(S) dt + d_s f. \quad (27) $$

**Remark 3.** ($d_s f$ is not directly related to the original dynamics). $d_s f$ is completely determined by $T_t^* \theta$ but is not directly related to the original dynamics Eq. (18). Therefore, once the expression of $T$ in Eq. (6) and the choice of $\theta$ are determined, the perturbation term $d_s f$ is prescribed. However, the choice of $\theta$ is up to the user and may then be related to the original dynamics.

**Remark 4.** In particular, there is no noise in the original dynamics (Eq. 18) which could be correlated with the noise of the resulting stochastic scheme (21). That is why the Itō lemma directly applies in the Taylor development (B4) of $f$ and then in the Eq. (22), leading to Eq. (23) and the final SPDE. Indeed, unlike the Itō–Wentzell formula (Kunita, 1997) – a cornerstone of the LU scheme – there is no additional cross-correlation term between $T_t^*$ and $\tilde{\theta}(t + \Delta t)$. The final SPDE (24) makes the link between the solution $\theta$ and the Brownian motions $\eta_i$ clear. But, at a given time step $t$,
since Eq. (18) has no noise term, $\tilde{\theta}(t + \Delta t)$ is correlated with the $t' \mapsto \eta_i(t')$ for $t' < t$ only and is independent of the new Brownian increment $\Delta \eta_i(t)$ generating $T_i$. Therefore, there is no cross-correlation term between $T_i^*$ and $\tilde{\theta}(t + \Delta t)$.

**Remark 5.** For a numerical implementation of our stochastic scheme, the time integration of the SPDE (25) may require a smaller time step $\Delta t$ than the time integration of the deterministic PDE (18) for two reasons. First, the available SDE (stochastic differential equation) time integration schemes are often less accurate than their deterministic counterparts. Secondly, the modified dynamics may involve additional Courant–Friedrichs–Lewy (CFL) constraints, related to for instance noise-induced diffusion.

**Example 3.2.1.** When $\theta = f$, Example (3.1.1),

$$ T_t^*\theta - \theta = (a^p \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f) \Delta t + e_i^p \partial_{x^p} f \Delta \eta_i. \quad (28) $$

This implies that

$$ d_s f = (a^p \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f) dt + e_i^p \partial_{x^p} f d\eta_i. \quad (29) $$

To physically interpret this equation, we rewrite

$$ \frac{d_s f}{dt} + V^p \partial_{x^p} f = \partial_{x^p} \left( \left( \frac{1}{2} e_i^p e_i^q \right) \partial_{x^q} f \right), \quad (30) $$

where

$$ V^p = -a^p + \frac{1}{2} \partial_{x^q} (e_i^p e_i^q) - e_i^p \frac{d\eta_i}{dt}. \quad (31) $$

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

If the original deterministic PDE (18) is an advection–diffusion equation, with advecting velocity $u$ and diffusion coefficient $D$, the final SPDE to simulate (Eq. 25) is now a stochastic advection–diffusion equation, with advecting velocity $u + V$ and diffusion matrix $DI_d + \frac{1}{2}e_i e_i^T$:

$$ \frac{df}{dt} + (u^p + V^p) \partial_{x^p} f = \partial_{x^p} \left( (D \delta_{pq} + \frac{1}{2} e_i^p e_i^q) \partial_{x^q} f \right). \quad (32) $$

This type of SPDE appears in the LU framework, detailed in Sect. 4.2.1.

**Example 3.2.2.** When $\theta = f dx^1 \wedge \dots \wedge dx^n$, Example 3.1.3,

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

This implies that

$$ d_s f = \left( \left( \partial_{x^p} a^p + \frac{1}{2} J_i f \right) f + \left( a^p + e_i^p \partial_{x^q} e_i^q \right) \partial_{x^p} f \right. \\ \left. + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f \right) dt \\ + \left( \partial_{x^p} e_i^p f + e_i^p \partial_{x^p} f \right) d\eta_i. \quad (34) $$

Rewritten, it leads to

$$ \frac{d_s f}{dt} + \partial_{x^p} (\tilde{V}^p f) = \partial_{x^p} \left( \left( \frac{1}{2} e_i^p e_i^q \right) \partial_{x^q} f \right), \quad (35) $$

where

$$ \tilde{V}^p = V^p - (e_i^p \partial_{x^q} e_i^q) = -a^p + \frac{1}{2} (\partial_{x^q} e_i^p e_i^q - e_i^p \partial_{x^q} e_i^q) - e_i^p \frac{d\eta_i}{dt}. \quad (36) $$

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

Recall, in fluid dynamics, that the Reynolds transport theorem provides an integral conservation equation for the transport of any conserved quantity within a fluid, connected to its corresponding differential equation. The Reynolds transport theorem is central to the LU setting. The present example thus already outlines a closed link between the proposed perturbation approach and the LU formulation. Accordingly, the SPDE (35) naturally appears in the LU framework, as detailed in Sect. 4.2.2.

**Example 3.2.3.** When $\theta = f^j dx^j$, Example 3.1.4,

$$ T_t^*\theta - \theta = \left\{ \left( a^p \partial_{x^p} f^j + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f^j + \partial_{x^j} a^p f^p \right. \right. \\ \left. \left. + \partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p \right) \Delta t \right. \\ \left. + (e_i^p \partial_{x^p} f^j + \partial_{x^j} e_i^p f^p) \Delta \eta_i \right\} dx^j. \quad (37) $$

For each $j$, the coefficients of $dx^j$ in $T_t^*\theta - \theta$ and those in $\theta$ can be compared, to lead to

$$ d_s f^j = (a^p \partial_{x^p} f^j + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f^j + \partial_{x^j} a^p f^p \\ + \partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p) dt \\ + (e_i^p \partial_{x^p} f^j + \partial_{x^j} e_i^p f^p) d\eta_i. \quad (38) $$

Regrouping the terms for physical interpretation, it reads

$$ \frac{d_s f^j}{dt} + V^p \partial_{x^p} f^j + \partial_{x^j} \left( -a^p - e_i^p \frac{d\eta_i}{dt} \right) f^p \\ - \partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p = \partial_{x^p} \left( \left( \frac{1}{2} e_i^p e_i^q \right) \partial_{x^q} f^j \right). \quad (39) $$
Two additional terms complete the advection–diffusion term. The first one, $\partial_{x^j} (-a^p - e_i^p \frac{d\eta_i}{dt}) f^p$, is reminiscent of the additional terms appearing in SALT momentum equations (Holm, 2015; Resseguier et al., 2020). The second term, $-\partial_{x^j} e_i^p e_i^q \partial_{x^q} f^p$, comes from the cross-correlation in Itô notation.

**Example 3.2.4.** When $\theta = f \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}$, Example 3.1.5,

$$ T_{t*}\theta - \theta = \left\{ \begin{aligned} & ((\partial_{x^p} a^p + \frac{1}{2} J_i) f + (-a^p + e_i^p \partial_{x^q} e_i^q) \\ & + \partial_{x^q} e_i^p e_i^q) \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f) \Delta t \\ & + (\partial_{x^p} e_i^p f - e_i^p \partial_{x^p} f) \Delta \eta_i \end{aligned} \right\} \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}. \quad (40) $$

This implies

$$ d_s f = \left( \begin{aligned} & (\partial_{x^p} a^p + \frac{1}{2} J_i) f + (-a^p + e_i^p \partial_{x^q} e_i^q) \\ & + \partial_{x^q} e_i^p e_i^q) \partial_{x^p} f + \frac{1}{2} e_i^p e_i^q \partial_{x^p} \partial_{x^q} f) dt \\ & + (\partial_{x^p} e_i^p f - e_i^p \partial_{x^p} f) d\eta_i. \end{aligned} \right) \quad (41) $$

It can then be verified that

$$ \frac{d_s f}{dt} + \partial_{x^p} \tilde{V}^p f - \tilde{V}^p \partial_{x^p} f = \partial_{x^p} \left( \left( \frac{1}{2} e_i^p e_i^q \right) \partial_{x^q} f \right), \quad (42) $$

where

$$ \tilde{V}^p = \tilde{V}^p - (e_i^p \partial_{x^q} e_i^q) = V^p - 2(e_i^p \partial_{x^q} e_i^q). \quad (43) $$

We recognize a diffusion term, $\partial_{x^p} (\frac{1}{2} e_i^p e_i^q) \partial_{x^q} f$, a velocity divergence term, $\partial_{x^p} \tilde{V}^p f$, and the advection term, $-\tilde{V}^p \partial_{x^p} f$. The divergence term is comparable to one appearing in the density equation.

However, the velocity fields appearing in the divergent and advecting terms do not coincide. Indeed, they are even opposite for divergence-free noise ($\partial_{x^q} e_i^q = 0$). This type of equation may appear uncommon but will be shown useful when applied to randomized thermal shallow water equations.

### 3.3 Conservation laws related to $d_s f$

A major advantage of the proposed perturbation scheme is to possibly prescribe $\theta$ to ensure that certain quantities are conserved. Define the discrete time version of $d_s f$ as

$$ \Delta_s f = \mathcal{M}^f(f) \Delta t + \mathcal{N}_i^f(f) \Delta \eta_i. \quad (44) $$

In general, conservation laws can be derived from the following two identities about the pull-back operator:

$$ (T_t^* \theta_1) \wedge (T_t^* \theta_2) = T_t^* (\theta_1 \wedge \theta_2) \quad (45) $$

$$ d T_t^* \theta = T_t d\theta, \quad (46) $$

where $d$ refers to the differential operator acting on differential forms. Hereafter, we present how to derive the conservation laws for two particular examples.

**Example 3.3.1.** Suppose $\theta_1 = f dx^1 \wedge \dots \wedge dx^n$ and define

$$ \hat{\theta}_1 = T_t^* \theta_1 \quad (47) $$

$$ \hat{f} = f + \Delta_s f. \quad (48) $$

Then $\hat{\theta}_1 = \hat{f} dx^1 \wedge \dots \wedge dx^n$. Therefore,

$$ \begin{aligned} \int_{\Omega} \hat{f} dx^1 \dots dx^n &= \int_{\Omega} \hat{\theta}_1 = \int_{\Omega} T_t^* \theta_1 = \int_{T_t(\Omega)} \theta_1 = \int_{\Omega} \theta_1 \\ &= \int_{\Omega} f dx^1 \dots dx^n. \end{aligned} \quad (49) $$

Equation (49) implies that the total integral of $f$ is not changed by the perturbation scheme. Next, suppose that $\theta_2 = g$ is a function. Similarly, we define

$$ \hat{\theta}_2 = T_t^* \theta_2 \quad (50) $$

$$ \hat{g} = g + \Delta_s g. \quad (51) $$

Applying Eq. (45),

$$ \begin{aligned} \int_{\Omega} \hat{f} \hat{g} dx^1 \dots dx^n &= \int_{\Omega} \hat{\theta}_1 \wedge \hat{\theta}_2 = \int_{\Omega} T_t^*(\theta_1 \wedge \theta_2) \\ &= \int_{T_t(\Omega)} \theta_1 \wedge \theta_2 = \int_{\Omega} \theta_1 \wedge \theta_2 = \int_{\Omega} fg dx^1 \dots dx^n. \end{aligned} \quad (52) $$

The total integral of $fg$ is thus also conserved by the perturbation scheme. Similarly for any integer $m \ge 0$, $fg^m$ is conserved by the perturbation scheme.

**Example 3.3.2.** Suppose $n = 2$ and $\theta = u dx + v dy$, where $u = (u, v)$ is the velocity field. The vorticity $\omega = \partial_x v - \partial_y u$ corresponds to the differential 2 form $d\theta$:

$$ d\theta = \omega dx^1 \wedge dx^2. \quad (53) $$

Define $\hat{\theta} := T_t^* \theta = \hat{u} dx^1 + \hat{v} dx^2$ and $\hat{\omega} = \partial_x \hat{v} - \partial_y \hat{u}$. Then $d\hat{\theta} = \hat{\omega} dx^1 \wedge dx^2$, and

$$ \begin{aligned} \int_{\Omega} \hat{\omega} dx^1 dx^2 &= \int_{\Omega} d\hat{\theta} = \int_{\Omega} d T_t^* \theta = \int_{\Omega} T_t^* d\theta \\ &= \int_{T_t(\Omega)} d\theta = \int_{\Omega} \omega dx^1 dx^2. \end{aligned} \quad (54) $$

Therefore, the vorticity is conserved by the perturbation scheme.

**Example 3.3.3.** Suppose $n = 3$ and $\theta = u dx + v dy + w dz$, where $u = (u, v, w)$ is the velocity field. The vorticity $\omega = (\partial_y w - \partial_z v, \partial_z u - \partial_x w, \partial_x v - \partial_y u)$ corresponds to the differential 2 form $d\theta$:

$$ \begin{aligned} d\theta &= (\partial_y w - \partial_z v) dy \wedge dz + (\partial_x v - \partial_y u) dz \wedge dx \\ &\quad + (\partial_x v - \partial_y u) dx \wedge dy. \end{aligned} \quad (55) $$
The helicity $Θ = u(\partial_y w - \partial_z v) + v(\partial_x v - \partial_y u) + w(\partial_x v - \partial_y u)$ corresponds to the differential 3 form:

$$d\theta \wedge \theta = (u(\partial_y w - \partial_z v) + v(\partial_x v - \partial_y u) + w(\partial_x v - \partial_y u))dx \wedge dy \wedge dz. \quad (56)$$

Similarly, we define $\hat{\Theta}$ by $d\hat{\theta} \wedge \hat{\theta} = \hat{\Theta} dx \wedge dy \wedge dz$. Then,

$$\begin{aligned} \int_{\Omega} \hat{\Theta} dx dy dz &= \int_{\Omega} d\hat{\theta} \wedge \hat{\theta} = \int_{\Omega} (dT_t^* \theta) \wedge (T_t^* \theta) \\ &= \int_{\Omega} (T_t^* d\theta) \wedge (T_t^* \theta) = \int_{\Omega} T_t^* (d\theta \wedge \theta) \\ &= \int_{T_t(\Omega)} d\theta \wedge \theta = \int_{\Omega} \Theta dx dy dz. \end{aligned} \quad (57)$$

Hence, in this case, the total amount of helicity is conserved.

**Example 3.3.4.** Suppose that $\theta_1 = f dx^1 \wedge \dots \wedge dx^n$ and that $\theta_2 = g \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}$. There exists a pairing $\langle \rangle$ for the differential $n$ forms and the contravariant $n$ vectors; i.e., $\langle \theta_1 \theta_2 \rangle = fg$ is a function on $\Omega$. Define

$$\hat{\theta}_1 = T_t^* \theta_1 = \hat{f} dx^1 \wedge \dots \wedge dx^n \quad (58)$$

$$\hat{\theta}_2 = (T_t^{-1})_* \theta_2 = \hat{g} \frac{\partial}{\partial x^1} \wedge \dots \wedge \frac{\partial}{\partial x^n}. \quad (59)$$

Then we have

$$\hat{f} \hat{g} (T_t^{-1}(x)) = \langle \hat{\theta}_1 \hat{\theta}_2 \rangle |_{T_t^{-1}(x)} = \langle \theta_1 \theta_2 \rangle |_x = fg(x), \quad (60)$$

and

$$\begin{aligned} \int_{\Omega} \hat{f}^2 \hat{g} dx^1 \dots dx^n &= \int_{\Omega} \langle \hat{\theta}_1 \hat{\theta}_2 \rangle \theta_1 = \int_{\Omega} \langle \theta_1 \theta_2 \rangle \theta_1 \\ &= \int_{\Omega} f^2 g dx^1 \dots dx^n. \end{aligned} \quad (61)$$

**Remark 6** (The conservation law of the perturbation scheme is independent of the conservation law of the original dynamical system). The derivation of Eqs. (49), (52), (54), (57), and (61) is based on the generic properties of the pull-back and push-forward operator of tensor fields. Since the choice of $\theta$ is not directly determined by the dynamical system, the conservation law of the perturbation scheme is independent of the original dynamical system. Recall that the perturbed forecast consists of two steps: Eqs. (20) and (21). The conservation law of the perturbation scheme implies that certain quantities are conserved in the second step. On the other hand, the original dynamical system (Eq. 20) might enjoy some other conservation law. If a quantity is conserved by both the original dynamical system and the perturbation scheme, then this quantity must be conserved by the final stochastic PDE. If a quantity is conserved by only one of Eqs. (20) and (21), then it can not be concluded that this quantity is conserved by the final SPDE.

## 4 Comparison with other perturbation schemes

In this section, we demonstrate that both the stochastic advection by Lie transport (SALT) equation (Holm, 2015) and the location uncertainty (LU) equation (Mémin, 2014; Resseguier et al., 2017a, 2020) can be recovered using the proposed perturbation scheme and properly choosing $\theta$ and the parameters $a, e_i$.

Note that the original LU paper (Mémin, 2014) assumed strong smoothness properties (finite variations in time) of the stochastic Navier–Stokes equations solution, to eventually remove the noises terms of this original Navier–Stokes equations under location uncertainty. Since Resseguier et al. (2017a, b), this assumption was removed, in order to keep the important noise terms. Accordingly, the original deterministic LU Navier–Stokes equations from Mémin (2014) have been referred to as pseudo-stochastic Navier–Stokes equations (Resseguier et al., 2021). Being deterministic, these pseudo-stochastic equations cannot be recovered by our stochastic scheme, whereas we can recover the stochastic LU Navier–Stokes equations that originated from Resseguier et al. (2017a).

### 4.1 Comparison with SALT equation

The original SALT equation (Holm, 2015) is derived based on a stochastically constrained variational principle $\delta S = 0$, for which

$$\left\{ \begin{array}{l} S(u, q) = \int \ell(u, q) dt \\ dq + \mathcal{L}_{dx_t} q = 0, \end{array} \right. \quad (62)$$

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

$$x_t(x) = x_0(x) + \int_0^t u(x, s) ds - \int_0^t e_i(x) \circ d\eta_i(s), \quad (63)$$

in which $u$ is the velocity vector field, and the $\circ$ means that the integral is defined in the Stratonovich sense, instead of in the Itô sense. Hence, $dx_t = u(x, t)dt - e_i \circ d\eta_i$ refers to an infinitesimal stochastic tangent field on the domain. Broadly speaking, we can express $dx_t = T_t(x) - x + udt$. Note the difference between Itô's notation and Stratonovich's notation; i.e., $e_i \circ d\eta_i \neq e_i d\eta_i$. Our expression of $T_t$ essentially follows Itô's notation, and $T_t(x) \neq x - e_i \Delta \eta_i$ in this subsection. Instead, it becomes $T_t(x) = x + \frac{1}{2} e_i^p \partial_{x_p} e_i \Delta t - e_i \Delta \eta_i$.

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

$$\mathcal{L}_{dx_t} q = d(i_{dx_t} q) + i_{dx_t} dq. \quad (64)$$
Essentially, the Lie derivative $\mathcal{L}_{dx_i} q$ corresponds to $T_t^* q - q + f^q(S)dt$, if we assume that the deterministic forecast of $q$ is simply the advection of $q$ by $u$. More generally, $\mathcal{L}_{dx_t - udt} q = T_t^* q - q$. Therefore, the SALT equation for $q$ is the same as our equation for $q$. We remark that Cartan's formula can not be directly applied to calculate the Lie derivative if the expression of $dx_t$ is in Itô's notation.

The SALT equation regarding the velocity $u$ comes from the first equation of Eq. (62). For most cases, the velocity $u$ is associated with the momentum, a differential 1 form $m = u^j dx^j = u^1 dx^1 + \dots + u^n dx^n$. In the examples discussed in Holm (2015), it is observed that, when the Lagrangian includes the kinetic energy, the stochastic noise contributes a term $\mathcal{L}_{dx_i} \theta$, where $\theta$ is a differential 1 form related to the momentum 1 form. For instance, $\theta = m$ in the example of “Stratonovich stochastic Euler-Poincaré flow” in Holm (2015), and $\theta = m + R^j dx^j$ in the example of “stochastic Euler-Boussinesq equations of a rotating stratified incompressible fluid” in Holm (2015). Already pointed out, the operator $\mathcal{L}_{dx_i}$ is closely related to $T_t^*$, and the momentum equation in SALT can be derived using our proposed scheme by properly choosing $\theta$.

Holm (2015) requires that $q$ is a differential form since Cartan's formula is only useful for differential forms $q$. This restriction can be relaxed by employing the original definition of Lie derivative with respect to a deterministic/stochastic flow of diffeomorphism discussed in Leon (2021), so that $\mathcal{L}_{dx_i} q$ can be generalized to the case where $q$ is a mixed-tensor field. This corresponds to Eq. (7) of the current paper.

Compared with Holm (2015) and Leon (2021), the proposed perturbation approach seems more flexible and does not have to rely on the Lagrangian mechanics. In particular, the velocity field can be associated with other tensor fields than the momentum 1 form. The perturbation, not directly related to the physics, can then be applied to any PDE. Moreover, our approach provides a new interpretation of $\mathcal{L}_{dx_t - udt}$ in terms of the optimal transportation associated with the infinitesimal forecast error at each time step. This interpretation certainly suggests practical numerical methods to infer $a, e_i$. Given a long sequence of reanalysis data or simulated high-resolution data, the one-step forecast can be evaluated using the low-resolution model, with the high-resolution state at each time step being the initial condition. $T_t$ is then estimated at each time step by comparing the low-resolution forecast and the high-resolution forecast. Finally, $a$ and $e_i$ could be learned from these samples of $T_t$.

## 4.2 Comparison with the LU equation

Mentioned above, the Reynolds transport theorem is central to the LU setting, and we already outlined a closed link between the proposed perturbation approach and the LU formulation. This link – related to differential $n$ forms – will be described precisely later in this subsection. But, before this,

we focus on another key ingredient of LU: the stochastic material derivative of functions (differential 0 forms).

### 4.2.1 0 forms in the LU framework

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

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

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

where $f$ can be any quantity that is assumed to be transported by the flow, i.e., $Df/Dt = 0$, where $D/Dt$ is the Itô material derivative. For instance, $f$ could be the velocity (dropping forces in the SPDE), the temperature, or the buoyancy. Compared to SALT notations, $-e_i d\eta_i$ is denoted $\boldsymbol{\sigma} d\boldsymbol{B} = \boldsymbol{\sigma}_{\bullet i} dB_i$. We refer to Resseguier et al. (2020, Appendix A) for the complete table of SALT-LU notation correspondences. Derived in Resseguier (2017, Appendix 10.1) and Resseguier et al. (2021, 6.1.3), we can rewrite it as

$$ \partial_t f + w_S \cdot \nabla f = \frac{1}{2}(\boldsymbol{\sigma}_{\bullet i} \cdot \nabla)(\boldsymbol{\sigma}_{\bullet i} \cdot \nabla f) - (\boldsymbol{\sigma} \dot{\boldsymbol{B}}) \cdot \nabla f, \quad (67) $$

$$ = -(\boldsymbol{\sigma} \circ \dot{\boldsymbol{B}}) \cdot \nabla f, \quad (68) $$

$$ w_S = w + w_S^c \quad (69) $$

$$ w_S^c = -\frac{1}{2}(\nabla \cdot \boldsymbol{a})^{\top} + \frac{1}{2}\boldsymbol{\sigma}(\nabla \cdot \boldsymbol{\sigma})^{\top}, \quad (70) $$

$$ = -\frac{1}{2}(\boldsymbol{\sigma}_{\bullet i} \cdot \nabla)\boldsymbol{\sigma}_{\bullet i}, \quad (71) $$

where $\boldsymbol{\sigma} \circ \dot{\boldsymbol{B}}$ is the Stratonovich noise of the SPDE, and $w$ and $w_S$ (denoted $u$ in the SALT framework) are the Itô drift and the Stratonovich drift of the fluid flow respectively. Separating the terms of the SPDE related to the deterministic dynamics from the term associated with the stochastic scheme results in

$$ d^{\text{LU}} f = g^f(S)dt + d_S^{\text{LU}} f, \quad (72) $$

where

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

$$ d_S^{\text{LU}} f = -w_S^c \cdot \nabla f dt + \frac{1}{2} (\boldsymbol{\sigma}_{\bullet i} \cdot \nabla)(\boldsymbol{\sigma}_{\bullet i} \cdot \nabla f) dt \\ - (\boldsymbol{\sigma} d\boldsymbol{B}) \cdot \nabla f. \quad (74) $$

Terms in Eqs. (65) and (66) translate to our notation in the following way:

$$ -w_S^c \cdot \nabla f dt = \frac{1}{2} e_i^q \partial_{x_q} e_i^p \partial_{x_p} f $$

$$ \frac{1}{2}(\boldsymbol{\sigma}_{\bullet i} \cdot \nabla)(\boldsymbol{\sigma}_{\bullet i} \cdot \nabla f) = \frac{1}{2} e_i^p \partial_{x_p} (e_i^q \partial_{x_q} f) $$

$$ = \frac{1}{2} (e_i^p \partial_{x_p} e_i^q \partial_{x_q} f + e_i^p e_i^q \partial_{x_p} \partial_{x_q} f) $$

$$ -\boldsymbol{\sigma} d\boldsymbol{B} \cdot \nabla f = e_i^p \partial_{x_p} f d\eta_i. $$
Hence,

$$d_s^{\text{LU}} f = (e_i^q \partial_{x_q} e_i^p \partial_{x_p} f + \frac{1}{2} e_i^p e_i^q \partial_{x_p} \partial_{x_q} f) \mathrm{d}t \\ + e_i^p \partial_{x_p} f \mathrm{d}\eta_i. \qquad (75)$$

Recall that Eq. (29) can be obtained by our perturbation scheme, while $f$ is associated with a differential 0 form. Direct calculation means that Eq. (75) coincides with Eq. (29) when

$$T_t(x) = x + e_i^q \partial_{x_q} e_i \Delta t + e_i \Delta \eta_i = x - w_S^c \Delta t \\ + (-w_S^c \Delta t - \sigma \Delta \mathbf{B}). \qquad (76)$$

The LU equation can thus be derived by choosing $\theta = f$ and $T_t$ by Eq. (76). At first glance, it does not seem straightforward to make such a choice. Nevertheless, it can be recognized that the term $(-w_S^c \Delta t - \sigma \Delta \mathbf{B}) = (\frac{1}{2} e_i^q \partial_{x_q} e_i \Delta t + e_i \Delta \eta_i)$ is the Itō noise and its Itō-to-Stratonovich correction. Hence, it corresponds to the Stratonovich noise $e_i \circ d\eta_i$ of the flow associated with $T_t$. The additional drift $-w_S^c \Delta t$ is different in nature. It is related to the advection correction $w_S^c \cdot \nabla f$ in the LU setting. Indeed, in the LU framework, the Itō drift, $w$, is seen as the resolved large-scale velocity. That is why, in this framework, the deterministic dynamics (74) involves the Itō drift, $w$. This is also the reason why, under the LU derivation, the advected velocity is assumed to be given by the Itō drift, $w$. It differs from the Stratonovich drift $w_S = w + w_S^c$, used as an advection velocity in the SALT approach or in Mikulevicius and Rozovskii (2004) (where the Stratonovich drift is denoted $u$). Interested readers are referred to Resseguier et al. (2020, Appendix A) for a discussion on these assumptions. Note however that in all these approaches, the advecting velocity is always the Stratonovich drift. This can be seen, e.g., in the Stratonovich form of LU Eq. (68).

To also understand Eq. (76), the inverse flow can be considered. According to Appendix (A),

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

Considering $T_t$ to represent how much the model forecast differs from the true forecast at every time step, $T_t^{-1}$ can be understood to represent how much the true forecast differs from the model forecast at each time step. Therefore, the LU equation can be derived using the proposed perturbation scheme, choosing $\theta = f$ and assuming that the true forecast differs from the model forecast by a displacement prescribed by Eq. (77).

### 4.2.2 $n$ forms in the LU framework

The LU physical justification relies on a stochastic interpretation of fundamental conservation laws, typically conservation of extensive properties (i.e., integrals of functions over a spatial volume) like momentum, mass, matter, and energy

(Resseguier et al., 2017a). These extensive properties can be expressed by integrals of differential $n$ forms. For instance, the mass and the momentum are integrals of the differential $n$ forms $\rho dx^1 \wedge \dots \wedge dx^n$ and $\rho w dx^1 \wedge \dots \wedge dx^n$, respectively. In the LU framework, a stochastic version of the Reynolds transport theorem (Resseguier et al., 2017a, Eq. 28) is used to deal with these differential $n$ forms $\theta = f dx^1 \wedge \dots \wedge dx^n$. Assuming an integral conservation $\frac{d}{dt} \int_{V(t)} f = 0$ on a spatial domain $V(t)$ transported by the flow, that theorem leads to the following SPDE:

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

where $D/Dt$ denotes the Itō material derivative. Here again, forcing terms are dropped for the sake of readability. This SPDE can be rewritten using the expression of that material derivative (Eqs. 9 and 10 of Resseguier et al., 2017a):

$$\partial_t f + \nabla \cdot (\mathbf{w}_S f) = \frac{1}{2} \nabla \cdot (a \nabla f) + \frac{1}{2} \nabla \cdot (\sigma_{\bullet i} (\nabla \cdot \sigma_{\bullet i})^T f) \\ - \nabla \cdot (\sigma \dot{\mathbf{B}} f) \qquad (79)$$

$$= \frac{1}{2} \nabla \cdot (\sigma_{\bullet i} (\nabla \cdot (\sigma_{\bullet i} f))^T) - \nabla \cdot (\sigma \dot{\mathbf{B}} f) \qquad (80)$$

$$= -\nabla \cdot (\sigma \circ \dot{\mathbf{B}} f). \qquad (81)$$

The original deterministic equation and stochastic perturbation correspond to

$$g^f(S) = -\nabla \cdot (\mathbf{w} f) \qquad (82)$$

$$d_s^{\text{LU}} f = (-\nabla \cdot (w_S^c f) + \frac{1}{2} \nabla \cdot (a \nabla f) + \frac{1}{2} \nabla \cdot (\sigma_{\bullet i} (\nabla \cdot \sigma_{\bullet i})^T f)) \mathrm{d}t - \nabla \cdot (\sigma \mathrm{d}\mathbf{B} f) \qquad (83)$$

$$= \nabla \cdot \left( \left( \left( \frac{1}{2} \nabla \cdot \mathbf{a} \right)^T \mathrm{d}t - \sigma \mathrm{d}\mathbf{B} \right) f \right) + \nabla \cdot \left( \frac{1}{2} \mathbf{a} \nabla f \right) \mathrm{d}t. \qquad (84)$$

Identifying $\mathbf{a} = \sigma_{\bullet i} \sigma_{\bullet i}^T = e_i e_i^T$ and $\sigma \dot{\mathbf{B}} = -e_i d\eta_i$, Eq. (35) corresponds to Example 3.2.2 about $n$ forms, with

$$\tilde{V} = -a^p + \frac{1}{2}(\partial_{x^q} e_i^p e_i^q - e_i^p \partial_{x^q} e_i^q) - e_i^p \frac{d\eta_i}{dt} \\ = -\left(\frac{1}{2} \nabla \cdot \mathbf{a}\right)^T + \sigma \dot{\mathbf{B}}; \qquad (85)$$

i.e.,

$$a^p = \partial_{x^q} (e_i^p e_i^q) - (e_i^p \partial_{x^q} e_i^q) = e_i^q \partial_{x^q} e_i^p. \qquad (86)$$

Again, the remapping is obtained,

$$T_t(x) = x + e_i^q \partial_{x_q} e_i \Delta t + e_i \Delta \eta_i = x - w_S^c \Delta t \\ + (-w_S^c \Delta t - \sigma \Delta \mathbf{B}), \qquad (87)$$
previously derived for the differential 0 form in LU framework (Eq. 76). Therefore, the proposed approach also generalizes the LU framework for n forms and its capacity – given by the Reynolds transport theorem – to deal with extensive properties.

**Remark 7.** For incompressible flows, the LU equation further imposes that

$$ \begin{cases} \nabla \cdot \boldsymbol{\sigma} = 0 \\ \nabla \cdot \nabla \cdot \boldsymbol{a} = 0. \end{cases} \qquad (88) $$

Translating it into our notation, it reads

$$ \begin{cases} \partial_{x_p} e_i^p = 0 & \text{for each } i \\ \partial_{x_p} \partial_{x_q} (e_i^p e_i^q) = 0. \end{cases} $$

Applying the result in Example 3.1.2, straightforward calculation means Eq. (88) is equivalent to $T_t^*\theta = \theta$ for $\theta = dx^1 \wedge \dots \wedge dx^n$. Such a result was expected since constraints (Eq. 88) are obtained from the LU density conservation.

## 5 A stochastic version of thermal shallow water equation

In this section, the proposed approach is applied to derive a stochastic version of thermal shallow water equation. Another stochastic version of thermal shallow water equation can be found in Holm and Luesink (2021). The thermal shallow water equation is derived in Warneford and Dellar (2013):

$$ \frac{\partial h}{\partial t} + \nabla \cdot (h \bar{u}) = 0, \qquad (89) $$

$$ \frac{\partial \Theta}{\partial t} + (\bar{u} \cdot \nabla) \Theta = -\kappa(h\Theta - h_0\Theta_0), \qquad (90) $$

$$ \frac{\partial \bar{u}}{\partial t} + (\bar{u} \cdot \nabla) \bar{u} + f \hat{z} \times \bar{u} = -\nabla(h\Theta) + \frac{1}{2} h \nabla \Theta. \qquad (91) $$

This model can be used to describe a two-layer system under equivalent barotropic approximation. The upper layer is active but with a spatiotemporal varying density $\rho(x, t)$, while the lower layer is quiescent with a fixed constant density $\rho_0$. The state variable $h$ represents the height of the active layer, and $\Theta = g(\rho_0 - \rho)/\rho_0$ is the density contrast. $\bar{u}$ is the averaged horizontal velocity of the active layer at each column. Note that $\rho < \rho_0$ (hence $\Theta > 0$) in the scenario of equivalent barotropic approximation (Warneford and Dellar, 2013).

Stated in Warneford and Dellar (2013), the following physical quantities are conserved up to the forcing.

$$ \text{Total energy: } E = \int_{\Omega} \frac{1}{2} (h |\bar{u}|^2 + h^2 \Theta) d^2 x. \qquad (92) $$

$$ \text{Total mass: } \mathcal{M} = \int_{\Omega} h d^2 x. \qquad (93) $$

$$ \text{Total momentum: } M = \int_{\Omega} h \bar{u} d^2 x. \qquad (94) $$

The objective is thus to choose proper tensor fields $\theta_{\bar{u}}$, $\theta_h$, and $\theta_{\Theta}$ for the state variables $\bar{u}$, $h$, and $\Theta$, respectively, so that $E, \mathcal{M}$, and $M$ are conserved by the perturbation scheme. Again, it must be emphasized that the conservation law of the perturbation scheme does not directly imply that the same quantities are conserved by the final SPDE.

The domain is 2-dimensional. To conserve mass, the only choice for $\theta_h$ is $\theta_h = hdx^1 \wedge dx^2$, which is a differential 2 form. It plays the role of density. In order to conserve the momentum, we need the momentum to be a differential 2 form as well. Hence we must choose $\theta_{\bar{u}}$ to be a function (differential 0 form). Therefore, the only choice for $\theta_{\bar{u}}$ is $\theta_{\bar{u}} = \bar{u}$. This choice of $\theta_{\bar{u}}$ and $\theta_h$ implies that $h|\bar{u}|^2$ also corresponds to a 2 form $|\bar{u}|^2 \theta_h$. Hence the kinetic energy is automatically conserved by the perturbation scheme. This means that if we want $E$ to be conserved, we must select $\theta_{\Theta}$ so that $h^2 \Theta$ corresponds to a differential 2 form. Note that $\theta_h$ is already a 2 form. We must thus select $\theta_{\Theta}$ so that $h \Theta$ corresponds to a function. The only choice for $\theta_{\Theta}$ is the contravariant tensor $\theta_{\Theta} = \Theta \frac{\partial}{\partial x^1} \wedge \frac{\partial}{\partial x^2}$. In this case, $h \Theta$ corresponds to the differential 0 form $\langle \theta_h \theta_{\Theta} \rangle = h \Theta$, where $\langle \cdot, \cdot \rangle$ in this section is the natural pairing of covariant n-tensor fields and contravariant n-tensor fields.

In sum, we have chosen the following tensor fields:

$$ \theta_h = h dx^1 \wedge dx^2 \qquad (95) $$

$$ \theta_{\bar{u}j} = \bar{u}^j \quad (\text{for } j = 1, 2) \qquad (96) $$

$$ \theta_{\Theta} = \Theta \frac{\partial}{\partial x^1} \wedge \frac{\partial}{\partial x^2}. \qquad (97) $$

For

$$ T_t(x) = x + a \Delta t + e_i \Delta \eta_i, \qquad (98) $$

we have

$$ T_t^{-1}(x) = x + (-a + e_i^p \partial_{x_p} e_i) \Delta t - e_i \Delta \eta_i. \qquad (99) $$

Then $T_t^*\theta_h$, $T_t^*\theta_{\bar{u}}$, and $(T_t^{-1})_*\theta_{\Theta}$ can be calculated following Examples 3.1.3, 3.1.1, and 3.1.5. This further implies $d_s h$, $d_s \bar{u}$, and $d_s \Theta$, as shown in Examples 3.2.2, 3.2.1, and 3.2.4. Note that $T_t^{-1}$ instead of $T_t$ is applied to $\theta_{\Theta}$, as shown in Eq. (7). Finally, we end up with the following SPDE:

$$ \begin{aligned} dh = & -\nabla(h\bar{u})dt + (h(\partial_{x_p}a^p + \frac{1}{2}J_i) + a^p \partial_{x_p}h \\ & + \frac{1}{2}e_i^p e_i^q \partial_{x_p} \partial_{x_q}h + \partial_{x_p}h e_i^p \partial_{x_q} e_i^q)dt \\ & + (h \partial_{x_p} e_i^p + \partial_{x_p} h e_i^p) d\eta_i \end{aligned} \qquad (100) $$
$$ 
\begin{aligned}
d\Theta = & -\{\bar{u} \cdot \nabla\}\Theta - \kappa(h\Theta - h_0\Theta_0)dt \\
& + (\Theta(-\partial_{x_p}a^p + \partial_{x_p}(\partial_{x_q}e_i e_i^q)^p + \frac{1}{2}J_i) \\
& + \partial_{x_p}\Theta a^p + \frac{1}{2}e_i^p e_i^q \partial_{x_p}\partial_{x_q}\Theta - \partial_{x_p}\Theta e_i^p \partial_{x_q}e_i^q)dt \\
& - (\Theta\partial_{x_p}e_i^p - \partial_{x_p}\Theta e_i^p)d\eta_i
\end{aligned} 
\quad (101) $$

$$ 
\begin{aligned}
d\bar{u}^j = & -\{\bar{u} \cdot \nabla\}\bar{u} - f\hat{z} \times \bar{u} - \nabla(h\Theta) + \frac{1}{2}h\nabla\Theta)^j dt \\
& + (\partial_{x_p}\bar{u}^j a^p + \frac{1}{2}e_i^p e_i^q \partial_{x_p}\partial_{x_q}\bar{u}^j)dt \\
& + \partial_{x_p}\bar{u}^j e_i^p d\eta_i,
\end{aligned} 
\quad (102) $$

where $J_i = \partial_{x_p} e_i^p \partial_{x_q} e_i^q - \partial_{x_q} e_i^p \partial_{x_p} e_i^q$. And the total mass, total momentum and the total energy shall all be conserved by the perturbation scheme.

# 6 Conclusions

The starting point of this work was to question how the location of the state variable can be consistently perturbed, motivated by Brenier's theorem (Brenier, 1991) which suggests that the difference of two density fields can be represented by a transport map $T$. Noting that optimal transportation has a clean representation in terms of differential $n$ forms, we proposed to perturb the "location" of the state variable $S$, at every forecast time step, by perturbing the corresponding differential $k$ forms $\theta$ by $\theta \leftarrow T_t^*\theta$, where $T_t$ is a random diffeomorphism which deviates from the identity map infinitesimally.

Under this framework, we end up with a stochastic PDE of the state variable $S$ in the form

$$ dS = f(S)dt + d_S S, \quad (103) $$

where $f(S)dt$ is the incremental of $S$ given by the original deterministic system. The term $d_S S$ is the additional stochastic incremental of $S$ caused by the perturbation scheme.

In this paper, we generalize this scheme to a mixed type of tensor fields $\theta$. A key point is indeed to link the state variable $S$ with some tensor field $\theta$. The choice of $\theta$ can then correspond to the conservation laws of certain quantities. We describe in detail how to calculate $T_t^*$ and $T_{t*}$ and present results for several examples corresponding to different choices of $\theta$. We also discussed about the conservation laws for these examples. We emphasize that Brenier's theorem merely serves as the motivation but not the theoretical foundation of the proposed scheme, since the optimality of the displacement vector field needs to be rigorously defined for general tensor fields $\theta$ that are not positive differential $n$ forms.

Interestingly, similarities and differences can be studied between the proposed perturbation scheme and the existing stochastic physical SALT and LU settings (Holm, 2015; Mémin, 2014; Resseguier et al., 2017a). In particular, both

SALT and LU equations can be recovered using a prescribed definition of the random diffeomorphism $T_t$ used by the perturbation scheme. For illustration, a stochastic version of the thermal shallow water equation is presented. Compared with SALT and LU settings (Holm, 2015; Mémin, 2014; Resseguier et al., 2017a), the proposed perturbation scheme does not directly rely on the physics. Hence it is more flexible and can be applied to any PDE. Yet, the proposed derivation also provides interesting means to interpret the operator $\mathcal{L}_{dx_t - udt}$, appearing in the SALT equation. In terms of the optimal transportation, this term represents the infinitesimal forecast error at every forecast time step.

In order to apply the proposed perturbation scheme to any specific model, the parameters $a$ and $e_i$ must be determined specifically. Hence it is necessary to learn these parameters from existing data, experimental runs, or additional physical considerations (Resseguier et al., 2020, 2021). We anticipate this framework naturally provides a new perspective on how to learn these parameters. Likely, this task will invoke the need of numerical algorithms to estimate the optimal transportation map for general differential $k$ forms or even mixed type of tensor fields. This will be the subject of future investigations.

# Appendix A: Calculation of $T_t^{-1}$

Suppose that

$$ T_t(x) = x + a\Delta t + e_i \Delta \eta_i. \quad (A1) $$

We assume that $T_t^{-1}$ has the following form of expression:

$$ T_t^{-1}(x) = x + z\Delta t + b_i \Delta \eta_i. \quad (A2) $$

Our goal is to find $z$ and $b_i$. Then we have

$$ 
\begin{aligned}
x = T_t(T_t^{-1}(x)) &= T_t(x + z\Delta t + b_i \Delta \eta_i) \\
&= x + z\Delta t + b_i \Delta \eta_i + a|_{x+z\Delta t+b_i \Delta \eta_i} \Delta t \\
&\quad + e_i|_{x+z\Delta t+b_i \Delta \eta_i} \Delta \eta_i.
\end{aligned} 
\quad (A3) $$

Similar to the derivation in Sect. 3.1, we apply Taylor expansion and Itô's lemma and drop the terms of the higher-order infinitesimal:

$$ 
\begin{aligned}
a|_{x+a\Delta t+b_i\Delta \eta_i}\Delta t &= a|_{x}\Delta t + o(\Delta t) \\
e_i|_{x+z\Delta t+b_i\Delta \eta_i}\Delta \eta_i &= e_i|_{x}\Delta \eta_i + e_{ip}b_i^p|_{x}\Delta t + o(\Delta t).
\end{aligned} 
\quad (A4) $$

Therefore,

$$ 
\begin{aligned}
x = T_t(T_t^{-1}(x)) &= x + (z + a + e_{ip}b_i^p)\Delta t \\
&\quad + (b_i + e_i)\Delta \eta_i + o(\Delta t).
\end{aligned} 
\quad (A5) $$

This implies that

$$ b_i + e_i = 0 \quad (A6) $$

$$ z + a + e_{ip}b_i^p = 0. \quad (A7) $$
Therefore,

$$b_i = -e_i \qquad (A8)$$

$$z = -a + e_i e_i^p, \qquad (A9)$$

or equivalently,

$$T_t^{-1}(x) = x + (-a + e_i e_i^p) \Delta t - e_i \Delta \eta_i. \qquad (A10)$$

## Appendix B: Derivation of $T_t^*\theta$

Given coordinates $(x^1, \dots, x^n)$, when $\theta$ is a differential $k$ form, it can be written as

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

Since $T_t^*$ is linear, we may assume that

$$\theta = f dx^{i_1} \wedge \dots \wedge dx^{i_k} \qquad (B2)$$

for some $1 \le i_1 < \dots < i_k \le n$. Let $T_t(x) = (T_t^1(x), \dots, T_t^n(x))$, then

$$(T_t^* \theta)(x) = f(T_t(x)) dT_t^{i_1} \wedge \dots \wedge dT_t^{i_k}. \qquad (B3)$$

We calculate $f(T_t(x))$ and $dT_t^{i_1} \wedge \dots \wedge dT_t^{i_k}$ separately. We denote $\Delta x = T_t(x) - x = a\Delta t + e_i \Delta \eta_i$ and $H_f$ the Hessian matrix of $f$. At a given time $t$, $f$ is assumed independent from the noise $\Delta \eta_i(t)$. Then

$$f(T_t(x)) = f(x + \Delta x) = f(x) + \langle \nabla f \Delta x \rangle + \frac{1}{2} (\Delta x)^\top H_f \Delta x + o((\Delta x)^2) \qquad (B4)$$

$$\begin{aligned} & = f(x) + \langle \nabla f \Delta t + e_i \Delta \eta_i \rangle + \frac{1}{2} e_i^\top H_f e_i (\Delta \eta_i)^2 \qquad (B5) \\ & + \mathcal{O}((\Delta t)^2) + \mathcal{O}(\Delta t \Delta \eta_i) + o((\Delta t)^2) + o((\Delta \eta_i)^2) \\ & \qquad + o(\Delta t \Delta \eta_i). \end{aligned}$$

According to Itô's lemma, $d\eta_i d\eta = dt$, and we can replace $(\Delta \eta_i)^2$ with $\Delta t$. Hence,

$$\begin{aligned} f(T_t(x)) = f(x) &+ \langle \nabla f a \rangle \Delta t + \langle \nabla f e_i \rangle \Delta \eta_i \\ &+ \frac{1}{2} e_i^\top H_f e_i \Delta t + o(\Delta t) \end{aligned} \qquad (B7)$$

$$\begin{aligned} & = f(x) + \left(\langle \nabla f a \rangle + \frac{1}{2} e_i H_f e_i \right) \Delta t + \langle \nabla f e_i \rangle \Delta \eta_i \\ & \quad + o(\Delta t). \end{aligned} \qquad (B8)$$

Next,

$$\begin{aligned} T_t^*(dx^{i_1} \wedge \dots \wedge dx^{i_k}) &= dT_t^{i_1} \wedge \dots \wedge dT_t^{i_k} \\ &= (dx^{i_1} + da^{i_1} \Delta t + de_i^{i_1} \Delta \eta_i) \wedge \dots \wedge \\ & \quad (dx^{i_k} + da^{i_k} \Delta t + de_i^{i_k} \Delta \eta_i). \end{aligned} \qquad (B9)$$

Note that $da^{ij}$ and $de_i^{ij}$ refer to the spatial differentiation. Again, we apply the “discrete version” of Itô’s rule $(\Delta \eta_i)^2 = \Delta t$ and collect all the terms of order $\mathcal{O}(\Delta t)$ and $\mathcal{O}(\Delta \eta_i)$:

$$\begin{aligned} T_t^*(dx^{i_1} \wedge \dots \wedge dx^{i_k}) = & dx^{i_1} \wedge \dots \wedge dx^{i_k} \\ & + (\sum_{s=1}^k dx^{i_1} \wedge \dots \wedge da^{is} \wedge \dots \wedge dx^{i_k}) \Delta t \\ & + (\sum_{s=1}^k dx^{i_1} \wedge \dots \wedge de_i^{is} \wedge \dots \wedge dx^{i_k}) \Delta \eta_i \\ & + (\sum_{s<r} dx^{i_1} \wedge \dots \wedge de_i^{is} \wedge \dots \wedge de_i^{ir}) \\ & \wedge \dots \wedge dx^{i_k}) \Delta t \\ & + o(\Delta t). \end{aligned} \qquad (B10)$$

According to the chain rule, $da^{is} = \partial_{x^j} a^{is} dx^j$, $de_i^{is} = \partial_{x^j} e_i^{is} dx^j$. Note that $\partial_{x^j} e_i^{is}$ refers to the $i$th component of $\partial_{x^j} e_i$, where $\partial_{x^j} e_i = \frac{\partial e_i}{\partial x^j}$ and $e_i(x) \in \mathbb{R}^n$ is the $i$th basis vector field of $T_t$. Hence,

$$\begin{aligned} T_t^*(dx^{i_1} \wedge \dots \wedge dx^{i_k}) &= dx^{i_1} \wedge \dots \wedge dx^{i_k} + (\sum_{s=1}^k \partial_{x^j} a^{is} dx^{i_1} \wedge \dots \wedge dx^j \\ & \qquad \wedge \dots \wedge dx^{i_k}) \Delta t \\ & + (\sum_{s=1}^k \partial_{x^j} e_i^{is} dx^{i_1} \wedge \dots \wedge dx^j \wedge \dots \wedge dx^{i_k}) \Delta \eta_i \\ & + (\sum_{s<r} \partial_{x^j} e_i^{is} \partial_{x^l} e_i^{ir} dx^{i_1} \wedge \dots \wedge dx^j \wedge \dots \wedge dx^l \\ & \qquad \wedge \dots \wedge dx^{i_k}) \Delta t \\ & + o(\Delta t). \end{aligned} \qquad (B11)$$
Combining Eqs. (B8) and (B11), with application of Itô's lemma, all terms of order $o(\Delta t)$ are then removed, to obtain

$$ 
\begin{aligned}
T_t^* \theta &= f(T_t(x)) T_t^*(dx^{i1} \wedge \dots \wedge dx^{ik}) \\
&= \theta + \{ (\nabla f a) + \frac{1}{2} e_i^\top H_f e_i \} dx^{i1} \wedge \dots \wedge dx^{in} \\
&\quad + \sum_{s=1}^k f \partial_{x^j} a^{is} dx^{i1} \wedge \dots dx^j \wedge \dots \wedge dx^{ik} \\
&\quad + \left( \sum_{s<r} f \partial_{x^j} e_i^{is} \partial_{x^l} e_i^{ir} dx^{i1} \wedge \dots \wedge dx^j \wedge \dots \wedge dx^l \wedge \dots \wedge dx^{ik} \right) \Delta t \\
&\quad + \left( \sum_{s=1}^k \{ (\nabla f e_i) \partial_{x^j} e_i^{is} dx^{i1} \wedge \dots \wedge dx^j \wedge \dots \wedge dx^{ik} \} \right) \Delta t \\
&\quad + \{ (\nabla f e_i) dx^{i1} \wedge \dots \wedge dx^{ik} \\
&\quad + \sum_{s=1}^k f \partial_{x^j} e_i^{is} dx^{i1} \\
&\quad \wedge \dots \wedge dx^j \wedge \dots \wedge dx^{ik} \} \Delta \eta_i \\
&\quad + o(\Delta t). 
\end{aligned} 
\quad (B12) 
$$

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

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

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

**Data availability.** No data sets were used in this article.

**Author contributions.** YZ contributed the main idea and mathematical derivation and wrote the first draft of the paper. VR contributed all the physical interpretations of the derived equations, a much more detailed introduction, and a more detailed comparison of the proposed scheme with the LU equation. BC contributed his insight and enthusiasm to this problem and worked on the writing of the draft. This work was done based on the fruitful discussion of all the three authors since a very early stage.

**Competing interests.** The contact author has declared that none of the authors has any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** The authors would like to express their gratitude towards Wei Pan, Darryl Holm, Dan Crisan, Long Li,

and Etienne Mémin for their patient explanation and insightful discussion. The research of Yicun Zhen was supported by the ANR Melody project when he was a postdoc at Ifremer. The research of Valentin Resseguier is supported by the company SCALIAN DS and by France Relance through the MORAANE project. The research of Bertrand Chapron is supported by ERC EU SYNERGY project no. 856408-STUOD and the ANR Melody project.

**Financial support.** This research has been supported by the Agence Nationale de la Recherche, Labex Immuno-Oncology (Melody), the European Research Council, H2020 (grant no. STUOD (856408)), the company SCALIAN DS and France Relance through the MORAANE project.

**Review statement.** This paper was edited by Jinqiao Duan and reviewed by Daniel Goodair and Baylor Fox-Kemper.

## References

Anderson, J. L.: An adaptive covariance inflation error correction algorithm for ensemble filters, Tellus A, 59, 210–224, 2007.

Berner, J., Shutts, G., Leutbecher, M., and Palmer, T.: A spectral stochastic kinetic energy backscatter scheme and its impact on flow-dependent predictability in the ECMWF ensemble prediction system, J. Atmos. Sci., 66, 603–626, 2009.

Brenier, Y.: Polar factorization and monotone rearrangement of vector-valued functions, Commun. Pur. Appl. Math., 44, 375–417, 1991.

Brzeźniak, Z., Capiński, M., and Flandoli, F.: Stochastic partial differential equations and turbulence, Math. Models Methods Appl. Sci., 1, 41–59, 1991.

Buizza, R., Miller, M., and Palmer, T.: Stochastic representation of model uncertainties in the ECMWF Ensemble Prediction System, Q. J. Roy. Meteor. Soc., 125, 2887–2908, 1999.

Chern, S. S., Chen, W. H., and Lam, K. S.: Lectures on Differential Geometry, World Scientific, https://doi.org/10.1142/3812, 1999.

Daum, F. E. and Huang, J.: Curse of dimensionality and particle filters, 2003 IEEE Aerospace Conference Proceedings (Cat. No.03TH8652), 4, 4_1979-4_1993, Big Sky, MT, USA, 8–15 March 2003, https://doi.org/10.1109/AERO.2003.1235126, 2003.

Flandoli, F.: The interaction between noise and transport mechanisms in PDEs, Milan J. Math., 79, 543–560, 2011.

Frank, J. E. and Gottwald, G. A.: Stochastic homogenization for an energy conserving multi-scale toy model of the atmosphere, Physica D, 254, 46–56, 2013.

Franzke, C., O’Kane, T., Berner, J., Williams, P., and Lucarini, V.: Stochastic climate theory and modeling, Wiley Interdiscip. Rev.: Clim. Change, 6, 63–78, https://doi.org/10.1002/wcc.318, 2015.

Gugole, F. and Franzke, C.: Numerical development and evaluation of an energy conserving conceptual stochastic climate model, Math. Clim. Weather Forecast., 5, 45–64, https://doi.org/10.1515/mcwf-2019-0004, 2019.

Harlim, J. and Majda, A. J.: Catastrophic filter divergence in filtering nonlinear dissipative systems, Commun. Math. Sci., 8, 27–43, 2010.
Hasselmann, K.: Stochastic climate models. Part I: theory, Tellus, 28, 473–485, 1976.

Holm, D. D.: Variational principles for stochastic fluid dynamics, Proc. Math. Phys. Eng. Sci., 471, https://doi.org/10.1098/rspa.2014.0963, 2015.

Holm, D. D. and Luesink, E.: Stochastic Wave–Current Interaction in Thermal Shallow Water Dynamics, J. Nonlinear Sci., 31, 29, https://doi.org/10.1007/s00332-021-09682-9, 2021.

Houtekamer, P. and Mitchell, H. L.: A Sequential Ensemble Kalman Filter for Atmospheric Data Assimilation, Mon. Wather Rev., 129, 123–137, 2001.

Jain, A., Timofeyev, I., and Vanden-Eijnden, E.: Stochastic mode-reduction in models with conservative fast sub-systems, arXiv [preprint], https://doi.org/10.48550/arXiv.1410.3004, 2014.

Kotsuki, S., Ota, Y., and Miyoshi, T.: Adaptive covariance relaxation methods for ensemble data assimilation: Experiments in the real atmosphere, Q. J. Roy. Meteor. Soc., 143, 2001–2015, 2017.

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

Leith, C.: Atmospheric predictability and two-dimensional turbulence, J. Atmos. Sci., 28, 145–161, 1971.

Leon, A. B. D.: On the effect of stochastic Lie transport noise on fluid dynamic equations, PhD thesis, Imperial College London, https://doi.org/10.25560/89936, 2021.

Leroux, S., Brankart, J.-M., Albert, A., Brodeau, L., Molines, J.-M., Jamet, Q., Le Sommer, J., Penduff, T., and Brasseur, P.: Ensemble quantification of short-term predictability of the ocean dynamics at a kilometric-scale resolution: a Western Mediterranean test case, Ocean Sci., 18, 1619–1644, https://doi.org/10.5194/os-18-1619-2022, 2022.

Leutbechner, M., Ollinaha, P., Lock, S.-J., Lang, S., Bechtold, P., Beljaars, A., Bozzo, A., Forbes, R., Haiden, T. Hogan, R., and Sandu, I.: Stochastic representations of model uncertainties in the IFS, in: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading, 2016.

Li, H., Kalnay, E., and Miyoshi, T.: Simultaneous estimation of covariance inflation and observation errors within an ensemble Kalman filter, Q. J. Roy. Meteor. Soc., 135, https://doi.org/10.1002/qj.371, 2009.

Majda, A. J., Timofeyev, I., and VandenEijnden, E.: Models for stochastic climate prediction, P. Natl. Acad. Sci. USA, 96 26, 14687–14691, 1999.

Mémin, E.: Fluid flow dynamics under location uncertainty, Geophys. Astro. Fluid, 108, 119–146, 2014.

Mikulevicius, R. and Rozovskii, B.: Stochastic Navier-Stokes Equations for Turbulent Flows, SIAM J. Math. Anal., 35, 1250–1310, 2004.

Miyoshi, T.: The Gaussian Approach to Adaptive Covariance Inflation and Its Implementation with the Local Ensemble Transform Kalman Filter, Mon. Wather Rev., 139, 1519–1535, 2011.

Orszag, S.: Analytical theories of turbulence, J. Fluid Mech., 41, 363–386, 1970.

Penland, C. and Sardeshmukh, P.: The optimal growth of tropical sea surface temperature anomalies, J. Climate, 8, 1999–2024, 1995.

Pope, S.: Lagrangian PDF methods for turbulent flows, Annu. Rev. Fluid Mech., 26, 23–63, 1994.

Poterjoy, J.: A Localized Particle Filter for High-Dimensional Non-linear Systems, Mon. Wather Rev., 144, 59–76, 2016.

Raanes, P. N., Bocquet, M., and Carrassi, A.: Adaptive covariance inflation in the ensemble Kalman filter by Gaussian scale mixtures, Q. J. Roy. Meteor. Soc., 145, 53–75, 2019.

Resseguier, V.: Mixing and fluid dynamics under location uncertainty, PhD thesis, Université Rennes 1, 2017.

Resseguier, V., Mémin, E., and Chapron, B.: Geophysical flows under location uncertainty, Part I Random transport and general models, Geophys. Astro. Fluid, 111, 149–176, 2017a.

Resseguier, V., Mémin, E., and Chapron, B.: Geophysical flows under location uncertainty, Part II Quasi-geostrophy and efficient ensemble spreading, Geophys. Astro. Fluid, 111, 177–208, 2017b.

Resseguier, V., Pan, W., and Fox-Kemper, B.: Data-driven versus self-similar parameterizations for stochastic advection by Lie transport and location uncertainty, Nonlin. Processes Geophys., 27, 209–234, https://doi.org/10.5194/npg-27-209-2020, 2020.

Resseguier, V., Li, L., Jouan, G., Dérian, P., Mémin, E., and Chapron, B.: New trends in ensemble forecast strategy: uncertainty quantification for coarse-grid computational fluid dynamics, Arch. Comput. Methods E., 28, 215–261, 2021.

Reynolds, C., Leutbecher, M., Batté, L., Chen, S., Christensen, H., Klasa, C., Pegion, P., Plant, B., Raynaud, L., Roberts, N., Sandu, I., Singleton, A., Sommer, M., Swinbank, R., Tennant, W., and Theis, S.: Reports from working group 3 : What are the pros/cons of existing model uncertainty schemes and how should these be measured?, in: ECMWF/WWRP Workshop: Model Uncertainty, ECMWF, Reading, 2016.

Sapsis, T. and Majda, A.: A statistically accurate modified quasilinear Gaussian closure for uncertainty quantification in turbulent dynamical systems, Physica D, 252, 34–45, 2013.

Schlee, F. H., Standish, C. J., and Toda, N. F.: Divergence in the Kalman Filter, AIAA J., 5, 1114–1120, 1966.

Tandeo, P., Ailliot, P., Bocquet, M., Carrassi, A., Carrassi, A., Miyoshi, T., Pulido, M., Pulido, M., and Zhen, Y.: A Review of Innovation-Based Methods to Jointly Estimate Model and Observation Error Covariance Matrices in Ensemble Data Assimilation, Mon. Wather Rev., 3973–3994, https://doi.org/10.1175/MWR-D-19-0240.1, 2020.

Tibshirani, R. and Knight, K.: The Covariance Inflation Criterion for Adaptive Model Selection, J. R. Stat. Soc. B, 61, 529–546, https://doi.org/10.1111/1467-9868.00191, 1999.

Warneford, E. S. and Dellar, P. J.: The quasi-geostrophic theory of the thermal shallow water equations, J. Fluid Mech., 723, 374–403, 2013.

Ying, Y. and Zhang, F.: An adaptive covariance relaxation method for ensemble data assimilation, Q. J. Roy. Meteor. Soc., 141, 2898–2906, https://doi.org/10.1002/qj.2576, 2015.

Zhen, Y. and Harlim, J.: Adaptive error covariances estimation methods for ensemble Kalman filters, J. Comput. Phys., 294, 619–638, 2015.