Physically constrained covariance inflation from location uncertainty

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 Majda, 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 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 (Penland 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 (ρ₁ and ρ₂) on a bounded domain Ω can be described by a transportation map: $T:\mathrm{\Omega }\to \mathrm{\Omega }$. More specifically, there exists a diffeomorphism T of Ω to transform ρ₁ to ρ₂ under the diffeomorphism T with a minimal cost. Broadly speaking, T can be interpreted as how much ρ₂ differs from ρ₁, and T operates as a location correction. Indeed, starting from the same initial condition ρ(t), suppose that ${\mathit{\rho }}_{\mathrm{1}}={\mathit{\rho }}^{\text{model}}(t+\mathrm{\Delta }t)$ is the model forecast and ${\mathit{\rho }}_{\mathrm{2}}=\mathit{\rho }(t+\mathrm{\Delta }t)$ is the true forecast. The additional uncertainty of ρ₁ 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=dimΩ. The statement “ρ₁ transforms to ρ₂ under the diffeomorphism T” is equivalent to the mathematical relation ${\mathit{\theta }}_{{\mathit{\rho }}_{\mathrm{1}}}=T^{*}{\mathit{\theta }}_{{\mathit{\rho }}_{\mathrm{2}}}$, where T^*, acting on all differential forms, is the pull-back operator induced by T or, equivalently, ${\mathit{\theta }}_{{\mathit{\rho }}_{\mathrm{2}}}=(T^{-\mathrm{1}}{)}^{*}{\mathit{\theta }}_{{\mathit{\rho }}_{\mathrm{1}}}$. Therefore, a random T (or equivalently, T⁻¹) 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:\mathrm{\Omega }\to \mathrm{\Omega }$.

• 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+Δt) based on the new value of S(t).

Associating S with different θ 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 θ 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 θ 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 ith 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 ∑_ie_i, and y_j refers to ∑_jy_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 Ω 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 \mathrm{0}$ and probability measures $\mathit{\mu },\mathit{\nu }\in \mathcal{P}\left(\mathrm{\Omega }\right)$,

$$\begin{array}{}\text{(1)}& \text{minimize }\mathbb{M}\left(T\right)={\int }_{\mathrm{\Omega }}c(x,T(x\left)\right)\mathrm{d}\mathit{\mu }\left(x\right)\end{array}$$

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

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

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\nu }=T_{\mathrm{no}.}\mathit{\mu }⟺g\left(x\right)=f\left(T\right(x\left)\right)\left|J_T\right(x\left)\right|,\end{array}$$

where J_T(x) refers to the Jacobian matrix of T at x. If we associate ν and μ with differential n forms ${\mathit{\theta }}_{\mathit{\nu }}=f\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\cdots }\wedge \mathrm{d}{x}^n$ and ${\mathit{\theta }}_{\mathit{\mu }}=g\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\cdots }\wedge \mathrm{d}{x}^n$, then

$$\begin{array}{}\text{(3)}& \mathit{\nu }=T_{\mathrm{no}.}\mathit{\mu }⟺{\mathit{\theta }}_{\mathit{\mu }}=T^{*}{\mathit{\theta }}_{\mathit{\nu }}.\end{array}$$

Brenier (1991) proved the existence and uniqueness of the solution to Monge's optimal transport problem for $c(x,y)=|x-y{|}^{\mathrm{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 μ and ν be measures with bounded smooth density on a bounded domain Ω⊂ℝⁿ. Let $c\mathrm{(}x\mathrm{,}y\mathrm{)}\mathrm{=}\mathrm{|}x\mathrm{-}y{\mathrm{|}}^{\mathrm{2}}$. Then there is a convex function $\mathit{\varphi }\mathrm{:}\mathrm{\Omega }\mathrm{\to }\mathbb{R}$, such that $\mathrm{(}\mathrm{\nabla }\mathit{\varphi }{\mathrm{)}}_{\mathrm{no}\mathrm{.}}\mathit{\mu }\mathrm{=}\mathit{\nu }$. And $\mathrm{\nabla }\mathit{\varphi }\mathrm{:}x\mathrm{\to }x\mathrm{+}\mathrm{\nabla }\mathit{\varphi }{\mathrm{|}}_x$, defined μ− almost everywhere, is the unique solution to Monge's optimal transport problem.

The convexity of ϕ implies that the map ∇ϕ 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 Ω. Let ρ denote the density field. Let ρ^model(t+Δt) and ρ^true(t+Δ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:\mathrm{\Omega }\to \mathrm{\Omega }$, so that

$$\begin{array}{}\text{(4)}& {\mathit{\rho }}^{\text{true}}(x,t+\mathrm{\Delta }t)={\mathit{\rho }}^{\text{model}}\left(T\right(x),t+\mathrm{\Delta }t)J_T\left(x\right).\end{array}$$

Note that the transportation T hereinafter is equivalent to the mapping T⁻¹ used in the Introduction. Equation (4) can further be written in terms of differential form. Let ${\mathit{\theta }}_{\mathit{\rho }}=\mathit{\rho }\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$, then Eq. (4) is equivalent to

$$\begin{array}{}\text{(5)}& T^{*}{\mathit{\theta }}_{\mathit{\rho }}^{\text{model}}(t+\mathrm{\Delta }t)={\mathit{\theta }}_{\mathit{\rho }}^{\text{true}}(t+\mathrm{\Delta }t).\end{array}$$

For general differential forms θ, 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:

$$\begin{array}{}\text{(6)}& T_t\left(x\right)=x+a(t,x)\mathrm{\Delta }t+e_i(t,x)\mathrm{\Delta }{\mathit{\eta }}_i\left(t\right),\end{array}$$

where $a(t,x),e_i(t,x)\in {\mathbb{R}}^n$, $\mathrm{\Delta }{\mathit{\eta }}_i\left(t\right)\sim \mathcal{N}(\mathrm{0},\mathrm{\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)Δt refers to a systematic location error, and e_iΔη_i refers to a random location error.

Stated in the Introduction, the state variable S must first be associated with a differential form θ. Then at every time step, T_t induces a perturbation of θ(t) by $\mathit{\theta }\left(t\right)\to T_t^{*}\mathit{\theta }\left(t\right)$. 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 $\mathit{\theta }=\mathit{\rho }\frac{\partial }{\partial x^{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \frac{\partial }{\partial x^n}$, where $\mathit{\{}\frac{\partial }{\partial x^i}{\mathit{\}}}_{i\le n}$ forms a global basis of the tangent field. Then T_t induces a perturbation of θ by $\mathit{\theta }\left(t\right)\to T_{t*}\mathit{\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 θ is a mixture of covariant and contravariant tensor fields, the perturbation scheme is slightly more complicated. Assume that $T_t:{\mathrm{\Omega }}_{\mathrm{1}}\to {\mathrm{\Omega }}_{\mathrm{2}}$ is a diffeomorphism, and $\mathit{\theta }=v\otimes \mathit{\omega }$, where v and ω are contravariant or covariant tensor fields respectively on Ω₂. Then $T_t^{*}\mathit{\omega }$ is a covariant tensor field on Ω₁. However, T_t can not directly induce a contravariant tensor field on Ω₁. In order to get a tensor field on Ω₁, we consider $T_t^{-\mathrm{1}}:{\mathrm{\Omega }}_{\mathrm{2}}\to {\mathrm{\Omega }}_{\mathrm{1}}$ and apply the push-forward operator on v. In sum, we may define the perturbation to be

$$\begin{array}{}\text{(7)}& \mathit{\theta }\left(t\right)\to \left(\right(T_t^{-\mathrm{1}}{)}_{*}v)\otimes (T_t^{*}\mathit{\omega }).\end{array}$$

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

3.1 Calculation of $T_t^{*}\mathit{\theta }$ (or T_t*θ)

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^{*}\mathit{\theta }$ (or T_t*θ), a Taylor expansion and It$\widehat{\text{o}}$'s lemma can be used.

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

$$\begin{array}{}\text{(8)}& \mathit{\theta }=\sum _{i_{\mathrm{1}}<\mathrm{\dots }<i_k}{f}^{i_{\mathrm{1}},\mathrm{\dots },i_k}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}.\end{array}$$

Then

$$\begin{array}{}\text{(9)}& T_t^{*}\mathit{\theta }=\sum _{i_{\mathrm{1}}<\mathrm{\dots }<i_k}{f}^{i_{\mathrm{1}},\mathrm{\dots },i_k}\left(T_t\right(x\left)\right)T_t^{*}(\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}).\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\displaystyle }{T}_t^{*}\mathit{\theta }=\mathit{\theta }+\mathcal{M}\left(\mathit{\theta }\right)\mathrm{\Delta }t+{\mathcal{N}}_i\left(\mathit{\theta }\right)\mathrm{\Delta }{\mathit{\eta }}_i,\end{array}$$

for some differential k forms ℳ(θ) and 𝒩_i(θ). Hereafter, several examples of $T_t^{*}\mathit{\theta }$ are presented.

The full derivation of these examples is skipped. We further express all the terms in coordinates. For instance, we replace 〈∇f〉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 $\mathit{\theta }=f\frac{\partial }{\partial x^{i_{\mathrm{1}}}}\wedge \mathrm{\dots }\wedge \frac{\partial }{\partial x^{i_k}}$ is a contravariant tensor field,

$$\begin{array}{}\text{(11)}& T_{t*}\mathit{\theta }=f\left(T_t^{-\mathrm{1}}\right(x\left)\right)T_{t*}\left({\displaystyle \frac{\partial }{\partial x^{i_{\mathrm{1}}}}}\wedge \mathrm{\dots }\wedge {\displaystyle \frac{\partial }{\partial x^{i_k}}}\right).\end{array}$$

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

Example 3.1.1 When θ=f is a function (differential 0 form),

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\left(T_t^{*}\mathit{\theta }\right)& =f+(a^j{\partial }_{x^j}f+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{\Delta }t\\ & +e_i^p{\partial }_{x^p}f\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}\end{array}$$

Example 3.1.2 When $\mathit{\theta }=\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{d}{x}^{\mathrm{2}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$,

$$\begin{array}{}\text{(13)}& T_t^{*}\mathit{\theta }=\mathit{\{}\mathrm{1}+({\partial }_{x^p}{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_i)\mathrm{\Delta }t+{\partial }_{x^p}{e}_i^p\mathrm{\Delta }{\mathit{\eta }}_i\mathit{\}}\mathit{\theta },\end{array}$$

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 $\mathit{\theta }=f\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{T}_t^{*}\mathit{\theta }=& \mathit{\{}f+(({\partial }_{x^p}{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_i)f+(a^p+e_i^p{\partial }_{x^q}{e}_i^q){\partial }_{x^p}f\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{\Delta }t\\ & +({\partial }_{x^p}{e}_i^{p}f+e_i^p{\partial }_{x^p}f)\mathrm{\Delta }{\mathit{\eta }}_i\mathit{\}}\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n.\end{array}\end{array}$$

Example 3.1.4 When θ=f^jdx^j (note that by the convention of notation, $f^j\mathrm{d}{x}^j={\sum }_{j=\mathrm{1}}^n{f}^j\mathrm{d}{x}^j$),

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{T}_t^{*}\mathit{\theta }=& \mathit{\{}{f}^j+(a^p{\partial }_{x^p}{f}^j+{\displaystyle \frac{\mathrm{1}}{\mathrm{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)\mathrm{\Delta }t\\ & +(e_i^p{\partial }_{x^p}{f}^j+{\partial }_{x^j}{e}_i^p{f}^p)\mathrm{\Delta }{\mathit{\eta }}_i\mathit{\}}\mathrm{d}{x}^j.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{T}_{t*}\mathit{\theta }=& \mathit{\{}f+(({\partial }_{x^p}{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{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+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{\Delta }t\\ & +({\partial }_{x^p}{e}_i^{p}f-e_i^p{\partial }_{x^p}f)\mathrm{\Delta }{\mathit{\eta }}_i\mathit{\}}{\displaystyle \frac{\partial }{\partial x^{\mathrm{1}}}}\wedge \mathrm{\dots }\wedge {\displaystyle \frac{\partial }{\partial x^n}}\end{array}.\end{array}$$

3.2 Derivation of the SPDE

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

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{\partial S}{\partial t}}=g\left(S\right).\end{array}$$

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

$$\begin{array}{}\text{(18)}& {\displaystyle }\mathrm{d}f=g^f\left(S\right)\mathrm{d}t.\end{array}$$

This implies a propagation equation for θ:

$$\begin{array}{}\text{(19)}& {\displaystyle }\mathrm{d}\mathit{\theta }=g^{\mathit{\theta }}\left(S\right)\mathrm{d}t.\end{array}$$

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

$$\begin{array}{}\text{(20)}& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t)=\mathit{\theta }\left(t\right)+g^{\mathit{\theta }}\left(S\right(t\left)\right)\mathrm{\Delta }t\text{(21)}& {\displaystyle }\mathit{\theta }(t+\mathrm{\Delta }t)=T_t^{*}\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t),\end{array}$$

with $T_t^{*}\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t)=\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t)+\mathcal{M}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right(t+\mathrm{\Delta }t\left)\right)\mathrm{\Delta }t+{\mathcal{N}}_i\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right(t+\mathrm{\Delta }t\left)\right)\mathrm{\Delta }{\mathit{\eta }}_i+o\left(\mathrm{\Delta }t\right)$ for some differential forms $\mathcal{M}\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)$ and ${\mathcal{N}}_i\left(\stackrel{\mathrm{̃}}{\mathit{\theta }}\right)$.

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

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{T}_t^{*}\stackrel{\mathrm{̃}}{\mathit{\theta }}& (t+\mathrm{\Delta }t)=\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t)+\left(\mathcal{M}\right(\mathit{\theta }\left(t\right))+\mathcal{O}(\mathrm{\Delta }t\left)\right)\mathrm{\Delta }t\\ & +\left({\mathcal{N}}_i\right(\mathit{\theta }\left(t\right))+\mathcal{O}(\mathrm{\Delta }t\left)\right)\mathrm{\Delta }{\mathit{\eta }}_i+o\left(\mathrm{\Delta }t\right)\\ =& \stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t)+\mathcal{M}\left(\mathit{\theta }\right(t\left)\right)\mathrm{\Delta }t+{\mathcal{N}}_i\left(\mathit{\theta }\right(t\left)\right)\mathrm{\Delta }{\mathit{\eta }}_i+o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

Therefore,

$$\begin{array}{}\text{(23)}& \begin{array}{rl}\mathit{\theta }(t+\mathrm{\Delta }t)& =\mathit{\theta }\left(t\right)+g^{\mathit{\theta }}\left(S\right(t\left)\right)\mathrm{\Delta }t+\mathcal{M}\left(\mathit{\theta }\right(t\left)\right)\mathrm{\Delta }t\\ & +{\mathcal{N}}_i\left(\mathit{\theta }\right(t\left)\right)\mathrm{\Delta }{\mathit{\eta }}_i+o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

This suggests the following stochastic propagation equation for θ:

$$\begin{array}{}\text{(24)}& \mathrm{d}\mathit{\theta }=g^{\mathit{\theta }}\left(S\right)\mathrm{d}t+\mathcal{M}\left(\mathit{\theta }\right)\mathrm{d}t+{\mathcal{N}}_i\left(\mathit{\theta }\right)\mathrm{d}{\mathit{\eta }}_i.\end{array}$$

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

$$\begin{array}{}\text{(25)}& \mathrm{d}f=g^f\left(S\right)\mathrm{d}t+{\mathcal{M}}^f\left(f\right)\mathrm{d}t+{\mathcal{N}}_i^f\left(f\right)\mathrm{d}{\mathit{\eta }}_i.\end{array}$$

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

$$\begin{array}{}\text{(26)}& {\mathrm{d}}_{s}f:={\mathcal{M}}^f\left(f\right)\mathrm{d}t+{\mathcal{N}}_i^f\left(f\right)\mathrm{d}{\mathit{\eta }}_i.\end{array}$$

Then Eq. (25) can be written as

$$\begin{array}{}\text{(27)}& \mathrm{d}f=g^f\left(S\right)\mathrm{d}t+{\mathrm{d}}_{s}f.\end{array}$$

Remark 3. (d_sf is not directly related to the original dynamics). d_sf is completely determined by $T_t^{*}\mathit{\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 θ are determined, the perturbation term d_sf is prescribed. However, the choice of θ 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 $\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t)$. The final SPDE (24) makes the link between the solution θ and the Brownian motions η_i clear. But, at a given time step t, since Eq. (18) has no noise term, $\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\Delta }t)$ is correlated with the $t^{\prime }\mapsto {\mathit{\eta }}_i\left(t^{\prime }\right)$ for $t^{\prime }<t$ only and is independent of the new Brownian increment Δη_i(t) generating T_t. Therefore, there is no cross-correlation term between $T_t^{*}$ and $\stackrel{\mathrm{̃}}{\mathit{\theta }}(t+\mathrm{\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 Δ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 θ=f, Example (3.1.1),

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{T}_t^{*}\mathit{\theta }-\mathit{\theta }& =(a^p{\partial }_{x^p}f+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{\Delta }t\\ & +e_i^p{\partial }_{x^p}f\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}\end{array}$$

This implies that

$$\begin{array}{}\text{(29)}& {\mathrm{d}}_{s}f=(a^p{\partial }_{x^p}f+{\displaystyle \frac{\mathrm{1}}{\mathrm{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}{\mathit{\eta }}_i.\end{array}$$

To physically interpret this equation, we rewrite

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{{\mathrm{d}}_{s}f}{\mathrm{d}t}}+V^p{\partial }_{x^p}f={\partial }_{x^p}\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q\right){\partial }_{x^q}f\right),\end{array}$$

where

$$\begin{array}{}\text{(31)}& V^p=-a^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\partial }_{x^q}\left(e_i^p{e}_i^q\right)-e_i^p{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_i}{\mathrm{d}t}}.\end{array}$$

Terms of advection and diffusion are recognized. The matrix $\frac{\mathrm{1}}{\mathrm{2}}{e}_i{e}_i^T$ is symmetric and non-negative and represents a diffusion matrix. The pth component of the advecting velocity V^p is composed of the drift −a^p, a correction $\frac{\mathrm{1}}{\mathrm{2}}{\partial }_{x^q}\left(e_i^p{e}_i^q\right)$, and a stochastic advecting velocity $-e_i^p\frac{\mathrm{d}{\mathit{\eta }}_i}{\mathrm{d}t}$.

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 $D{I}_d+\frac{\mathrm{1}}{\mathrm{2}}{e}_i{e}_i^T$:

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\mathrm{d}f}{\mathrm{d}t}}+(u^p+V^p){\partial }_{x^p}f={\partial }_{x^p}\left((D{\mathit{\delta }}_{pq}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q){\partial }_{x^q}f\right).\end{array}$$

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

Example 3.2.2. When $\mathit{\theta }=f\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$, Example 3.1.3,

$$\begin{array}{}\text{(33)}& \begin{array}{rl}{T}_t^{*}& \mathit{\theta }-\mathit{\theta }=\mathit{\left\{}\right(({\partial }_{x^p}{a}^p\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_i)f+(a^p+e_i^p{\partial }_{x^q}{e}_i^q){\partial }_{x^p}f+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{\Delta }\\ & +({\partial }_{x^p}{e}_i^{p}f+e_i^p{\partial }_{x^p}f)\mathrm{\Delta }{\mathit{\eta }}_i\mathit{\}}\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n.\end{array}\end{array}$$

This implies that

$$\begin{array}{}\text{(34)}& \begin{array}{rl}{\mathrm{d}}_{s}f=& \left(\right({\partial }_{x^p}{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_i)f+(a^p+e_i^p{\partial }_{x^q}{e}_i^q){\partial }_{x^p}f\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{d}t\\ & +({\partial }_{x^p}{e}_i^{p}f+e_i^p{\partial }_{x^p}f)\mathrm{d}{\mathit{\eta }}_i.\end{array}\end{array}$$

Rewritten, it leads to

$$\begin{array}{}\text{(35)}& {\displaystyle \frac{{\mathrm{d}}_{s}f}{\mathrm{d}t}}+{\partial }_{x^p}\left({\stackrel{\mathrm{̃}}{V}}^{p}f\right)={\partial }_{x^p}\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q\right){\partial }_{x^q}f\right),\end{array}$$

where

$$\begin{array}{}\text{(36)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{V}}^p& =V^p-\left(e_i^p{\partial }_{x^q}{e}_i^q\right)=-a^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\partial }_{x^q}{e}_i^p{e}_i^q-e_i^p{\partial }_{x^q}{e}_i^q)\\ & -e_i^p{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_i}{\mathrm{d}t}}.\end{array}\end{array}$$

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_t(x)) and $T_t^{*}(\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^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 θ=f^jdx^j, Example 3.1.4,

$$\begin{array}{}\text{(37)}& \begin{array}{rl}{T}_t^{*}\mathit{\theta }-\mathit{\theta }=& \mathit{\left\{}\right(a^p{\partial }_{x^p}{f}^j+{\displaystyle \frac{\mathrm{1}}{\mathrm{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)\mathrm{\Delta }t\\ & +(e_i^p{\partial }_{x^p}{f}^j+{\partial }_{x^j}{e}_i^p{f}^p)\mathrm{\Delta }{\mathit{\eta }}_i\mathit{\}}\mathrm{d}{x}^j.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{\mathrm{d}}_s{f}^j=& (a^p{\partial }_{x^p}{f}^j+{\displaystyle \frac{\mathrm{1}}{\mathrm{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)\mathrm{d}t\\ & +(e_i^p{\partial }_{x^p}{f}^j+{\partial }_{x^j}{e}_i^p{f}^p)\mathrm{d}{\mathit{\eta }}_i.\end{array}\end{array}$$

Regrouping the terms for physical interpretation, it reads

$$\begin{array}{}\text{(39)}& \begin{array}{rl}{\displaystyle \frac{{\mathrm{d}}_s{f}^j}{\mathrm{d}t}}& +V^p{\partial }_{x^p}{f}^j+{\partial }_{x^j}\left(-a^p-e_i^p{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_i}{\mathrm{d}t}}\right)f^p\\ & -{\partial }_{x^j}{e}_i^p{e}_i^q{\partial }_{x^q}{f}^p={\partial }_{x^p}\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q\right){\partial }_{x^q}{f}^j\right).\end{array}\end{array}$$

Two additional terms complete the advection–diffusion term. The first one, ${\partial }_{x^j}\left(-a^p-e_i^p\frac{\mathrm{d}{\mathit{\eta }}_i}{\mathrm{d}t}\right)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 $\mathit{\theta }=f\frac{\partial }{\partial x^{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \frac{\partial }{\partial x^n}$, Example 3.1.5,

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{T}_{t*}\mathit{\theta }-& \mathit{\theta }=\mathit{\left\{}\right(({\partial }_{x^p}{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{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+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{\Delta }t\\ & +({\partial }_{x^p}{e}_i^{p}f-e_i^p{\partial }_{x^p}f)\mathrm{\Delta }{\mathit{\eta }}_i\mathit{\}}{\displaystyle \frac{\partial }{\partial x^{\mathrm{1}}}}\wedge \mathrm{\dots }\wedge {\displaystyle \frac{\partial }{\partial x^n}}.\end{array}\end{array}$$

This implies

$$\begin{array}{}\text{(41)}& \begin{array}{rl}{\mathrm{d}}_{s}f=& \left(\right({\partial }_{x^p}{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{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+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x^p}{\partial }_{x^q}f)\mathrm{d}t\\ & +({\partial }_{x^p}{e}_i^{p}f-e_i^p{\partial }_{x^p}f)\mathrm{d}{\mathit{\eta }}_i.\end{array}\end{array}$$

It can then be verified that

$$\begin{array}{}\text{(42)}& {\displaystyle \frac{{\mathrm{d}}_{s}f}{\mathrm{d}t}}+{\partial }_{x^p}{\stackrel{\mathrm{̃}}{V}}^{p}f-{\stackrel{\mathrm{̃}}{\stackrel{\mathrm{̃}}{V}}}^p{\partial }_{x^p}f={\partial }_{x^p}\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q\right){\partial }_{x^q}f\right),\end{array}$$

where

$$\begin{array}{}\text{(43)}& {\stackrel{\mathrm{̃}}{\stackrel{\mathrm{̃}}{V}}}^p={\stackrel{\mathrm{̃}}{V}}^p-\left(e_i^p{\partial }_{x^q}{e}_i^q\right)=V^p-\mathrm{2}\left(e_i^p{\partial }_{x^q}{e}_i^q\right).\end{array}$$

We recognize a diffusion term, ${\partial }_{x^p}\left(\left(\frac{\mathrm{1}}{\mathrm{2}}{e}_i^p{e}_i^q\right){\partial }_{x^q}f\right)$, a velocity divergence term, ${\partial }_{x^p}{\stackrel{\mathrm{̃}}{V}}^{p}f$, and the advection term, $-{\stackrel{\mathrm{̃}}{\stackrel{\mathrm{̃}}{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=\mathrm{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_sf

A major advantage of the proposed perturbation scheme is to possibly prescribe θ to ensure that certain quantities are conserved. Define the discrete time version of d_sf as

$$\begin{array}{}\text{(44)}& {\mathrm{\Delta }}_{s}f={\mathcal{M}}^f\left(f\right)\mathrm{\Delta }t+{\mathcal{N}}_i^f\left(f\right)\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}$$

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

$$\begin{array}{}\text{(45)}& {\displaystyle }\left(T_t^{*}{\mathit{\theta }}_{\mathrm{1}}\right)\wedge \left(T_t^{*}{\mathit{\theta }}_{\mathrm{2}}\right)=T_t^{*}({\mathit{\theta }}_{\mathrm{1}}\wedge {\mathit{\theta }}_{\mathrm{2}})\text{(46)}& {\displaystyle }\mathrm{d}{T}_t^{*}\mathit{\theta }=T_{t}d\mathit{\theta },\end{array}$$

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 ${\mathit{\theta }}_{\mathrm{1}}=f\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$ and define

$$\begin{array}{}\text{(47)}& {\displaystyle }{\widehat{\mathit{\theta }}}_{\mathrm{1}}=T_t^{*}{\mathit{\theta }}_{\mathrm{1}}\text{(48)}& {\displaystyle }\widehat{f}=f+{\mathrm{\Delta }}_{s}f.\end{array}$$

Then ${\widehat{\mathit{\theta }}}_{\mathrm{1}}=\widehat{f}\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$. Therefore,

$$\begin{array}{}\text{(49)}& \begin{array}{rl}& \underset{\mathrm{\Omega }}{\int }\widehat{f}\mathrm{d}{x}^{\mathrm{1}}\mathrm{\dots }\mathrm{d}{x}^n=\underset{\mathrm{\Omega }}{\int }{\widehat{\mathit{\theta }}}_{\mathrm{1}}=\underset{\mathrm{\Omega }}{\int }{T}_t^{*}{\mathit{\theta }}_{\mathrm{1}}=\underset{T_t\left(\mathrm{\Omega }\right)}{\int }{\mathit{\theta }}_{\mathrm{1}}=\underset{\mathrm{\Omega }}{\int }{\mathit{\theta }}_{\mathrm{1}}\\ =& \underset{\mathrm{\Omega }}{\int }f\mathrm{d}{x}^{\mathrm{1}}\mathrm{\dots }\mathrm{d}{x}^n.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(50)}& {\displaystyle }{\widehat{\mathit{\theta }}}_{\mathrm{2}}=T_t^{*}{\mathit{\theta }}_{\mathrm{2}}\text{(51)}& {\displaystyle }\widehat{g}=g+{\mathrm{\Delta }}_{s}g.\end{array}$$

Applying Eq. (45),

$$\begin{array}{}\text{(52)}& \begin{array}{rl}\underset{\mathrm{\Omega }}{\int }& \widehat{f}\widehat{g}\mathrm{d}{x}^{\mathrm{1}}\mathrm{\dots }\mathrm{d}{x}^n=\underset{\mathrm{\Omega }}{\int }{\widehat{\mathit{\theta }}}_{\mathrm{1}}\wedge {\widehat{\mathit{\theta }}}_{\mathrm{2}}=\underset{\mathrm{\Omega }}{\int }{T}_t^{*}({\mathit{\theta }}_{\mathrm{1}}\wedge {\mathit{\theta }}_{\mathrm{2}})\\ & =\underset{T_t\left(\mathrm{\Omega }\right)}{\int }{\mathit{\theta }}_{\mathrm{1}}\wedge {\mathit{\theta }}_{\mathrm{2}}=\underset{\mathrm{\Omega }}{\int }{\mathit{\theta }}_{\mathrm{1}}\wedge {\mathit{\theta }}_{\mathrm{2}}=\underset{\mathrm{\Omega }}{\int }fg\mathrm{d}{x}^{\mathrm{1}}\mathrm{\dots }.\mathrm{d}{x}^n.\end{array}\end{array}$$

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

Example 3.3.2. Suppose n=2 and $\mathit{\theta }=u\mathrm{d}x+v\mathrm{d}y$, where $\mathit{u}=(u,v)$ is the velocity field. The vorticity $\mathit{\omega }={\partial }_{x}v-{\partial }_{y}u$ corresponds to the differential 2 form dθ:

$$\begin{array}{}\text{(53)}& \mathrm{d}\mathit{\theta }=\mathit{\omega }\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{d}{x}^{\mathrm{2}}.\end{array}$$

Define $\widehat{\mathit{\theta }}:=T_t^{*}\mathit{\theta }=\widehat{u}\mathrm{d}{x}^{\mathrm{1}}+\widehat{v}\mathrm{d}{x}^{\mathrm{2}}$ and $\widehat{\mathit{\omega }}={\partial }_x\widehat{v}-{\partial }_y\widehat{u}$. Then $\mathrm{d}\widehat{\mathit{\theta }}=\widehat{\mathit{\omega }}\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{d}{x}^{\mathrm{2}}$, and

$$\begin{array}{}\text{(54)}& \begin{array}{rl}\underset{\mathrm{\Omega }}{\int }& \widehat{\mathit{\omega }}\mathrm{d}{x}^{\mathrm{1}}\mathrm{d}{x}^{\mathrm{2}}=\underset{\mathrm{\Omega }}{\int }d\widehat{\mathit{\theta }}=\underset{\mathit{\omega }}{\int }\mathrm{d}{T}_t^{*}\mathit{\theta }=\underset{\mathrm{\Omega }}{\int }{T}_t^{*}\mathrm{d}\mathit{\theta }\\ & =\underset{T_t\left(\mathrm{\Omega }\right)}{\int }\mathrm{d}\mathit{\theta }=\underset{\mathrm{\Omega }}{\int }\mathit{\omega }\mathrm{d}{x}^{\mathrm{1}}\mathrm{d}{x}^{\mathrm{2}}.\end{array}\end{array}$$

Therefore, the vorticity is conserved by the perturbation scheme.

Example 3.3.3. Suppose n=3 and $\mathit{\theta }=u\mathrm{d}x+v\mathrm{d}y+w\mathrm{d}z$, where $\mathit{u}=(u,v,w)$ is the velocity field. The vorticity $\mathit{\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θ:

$$\begin{array}{}\text{(55)}& \begin{array}{rl}\mathrm{d}\mathit{\theta }& =({\partial }_{y}w-{\partial }_{z}v)\mathrm{d}y\wedge \mathrm{d}z+({\partial }_{x}v-{\partial }_{y}u)\mathrm{d}z\wedge \mathrm{d}x\\ & +({\partial }_{x}v-{\partial }_{y}u)\mathrm{d}x\wedge \mathrm{d}y.\end{array}\end{array}$$

The helicity $\mathrm{\Theta }=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:

$$\begin{array}{}\text{(56)}& \begin{array}{rl}\mathrm{d}\mathit{\theta }\wedge \mathit{\theta }& =\left(u\right({\partial }_{y}w-{\partial }_{z}v)+v({\partial }_{x}v-{\partial }_{y}u)\\ & +w({\partial }_{x}v-{\partial }_{y}u))\mathrm{d}x\wedge \mathrm{d}y\wedge \mathrm{d}z.\end{array}\end{array}$$

Similarly, we define $\widehat{\mathrm{\Theta }}$ by $\mathrm{d}\widehat{\mathit{\theta }}\wedge \widehat{\mathit{\theta }}=\widehat{\mathrm{\Theta }}\mathrm{d}x\wedge \mathrm{d}y\wedge \mathrm{d}z$. Then,

$$\begin{array}{}\text{(57)}& \begin{array}{rl}\underset{\mathrm{\Omega }}{\int }\widehat{\mathrm{\Theta }}\mathrm{d}x\mathrm{d}y\mathrm{d}z& =\underset{\mathrm{\Omega }}{\int }\mathrm{d}\widehat{\mathit{\theta }}\wedge \widehat{\mathit{\theta }}=\underset{\mathrm{\Omega }}{\int }\left(\mathrm{d}{T}_t^{*}\mathit{\theta }\right)\wedge \left(T_t^{*}\mathit{\theta }\right)\\ =& \underset{\mathrm{\Omega }}{\int }\left(T_t^{*}\mathrm{d}\mathit{\theta }\right)\wedge \left(T_t^{*}\mathit{\theta }\right)=\underset{\mathrm{\Omega }}{\int }{T}_t^{*}(\mathrm{d}\mathit{\theta }\wedge \mathit{\theta })\\ & =\underset{T_t\left(\mathrm{\Omega }\right)}{\int }\mathrm{d}\mathit{\theta }\wedge \mathit{\theta }=\underset{\mathrm{\Omega }}{\int }\mathrm{\Theta }\mathrm{d}x\mathrm{d}y\mathrm{d}z.\end{array}\end{array}$$

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

Example 3.3.4. Suppose that ${\mathit{\theta }}_{\mathrm{1}}=f\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$ and that ${\mathit{\theta }}_{\mathrm{2}}=g\frac{\partial }{\partial x^{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \frac{\partial }{\partial x^n}$. There exists a pairing 〈〉 for the differential n forms and the contravariant n vectors; i.e., 〈θ₁θ₂〉=fg is a function on Ω. Define

$$\begin{array}{}\text{(58)}& {\displaystyle }{\widehat{\mathit{\theta }}}_{\mathrm{1}}=T_t^{*}{\mathit{\theta }}_{\mathrm{1}}=\widehat{f}\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n\text{(59)}& {\displaystyle }{\widehat{\mathit{\theta }}}_{\mathrm{2}}=(T_t^{-\mathrm{1}}{)}_{*}{\mathit{\theta }}_{\mathrm{2}}=\widehat{g}{\displaystyle \frac{\partial }{\partial x^{\mathrm{1}}}}\wedge \mathrm{\dots }\wedge {\displaystyle \frac{\partial }{\partial x^n}}.\end{array}$$

Then we have

$$\begin{array}{}\text{(60)}& \widehat{f}\widehat{g}\left(T_t^{-\mathrm{1}}\right(x\left)\right)=\langle {\widehat{\mathit{\theta }}}_{\mathrm{1}}{\widehat{\mathit{\theta }}}_{\mathrm{2}}\rangle {|}_{T_t^{-\mathrm{1}}\left(x\right)}=\langle {\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\rangle {|}_x=fg\left(x\right),\end{array}$$

and

$$\begin{array}{}\text{(61)}& \begin{array}{rl}\underset{\mathrm{\Omega }}{\int }{\widehat{f}}^{\mathrm{2}}\widehat{g}\mathrm{d}{x}^{\mathrm{1}}\mathrm{\dots }\mathrm{d}{x}^n& =\underset{\mathrm{\Omega }}{\int }\langle {\widehat{\mathit{\theta }}}_{\mathrm{1}}{\widehat{\mathit{\theta }}}_{\mathrm{2}}\rangle {\mathit{\theta }}_{\mathrm{1}}=\underset{\mathrm{\Omega }}{\int }\langle {\mathit{\theta }}_{\mathrm{1}}{\mathit{\theta }}_{\mathrm{2}}\rangle {\mathit{\theta }}_{\mathrm{1}}\\ & =\underset{\mathrm{\Omega }}{\int }{f}^{\mathrm{2}}g\mathrm{d}{x}^{\mathrm{1}}\mathrm{\dots }\mathrm{d}{x}^n.\end{array}\end{array}$$

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 θ 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 θ 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 δS=0, for which

$$\begin{array}{}\text{(62)}& \left\{\begin{array}{ll}& S(u,q)=\int \mathrm{\ell }(u,q)\mathrm{d}t\\ & \mathrm{d}q+{\mathfrak{L}}_{\mathrm{d}{x}_t}q=\mathrm{0},\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(63)}& x_t\left(x\right)=x_{\mathrm{0}}\left(x\right)+\underset{\mathrm{0}}{\overset{t}{\int }}u(x,s)\mathrm{d}s-\underset{\mathrm{0}}{\overset{t}{\int }}{e}_i\left(x\right)\circ \mathrm{d}{\mathit{\eta }}_i\left(s\right),\end{array}$$

in which u is the velocity vector field, and the ∘ means that the integral is defined in the Stratonovich sense, instead of in the Itô sense. Hence, $\mathrm{d}{x}_t=u(x,t)\mathrm{d}t-e_i\circ \mathrm{d}{\mathit{\eta }}_i$ refers to an infinitesimal stochastic tangent field on the domain. Broadly speaking, we can express $\mathrm{d}{x}_t=T_t\left(x\right)-x+u\mathrm{d}t$. Note the difference between Itô's notation and Stratonovich's notation; i.e., $e_i\circ \mathrm{d}{\mathit{\eta }}_i\ne e_i\mathrm{d}{\mathit{\eta }}_i$. Our expression of T_t essentially follows Itô's notation, and $T_t\left(x\right)\ne x-e_i\mathrm{\Delta }{\mathit{\eta }}_i$ in this subsection. Instead, it becomes $T_t\left(x\right)=x+\frac{\mathrm{1}}{\mathrm{2}}{e}_i^p{\partial }_{x_p}{e}_i\mathrm{\Delta }t-e_i\mathrm{\Delta }{\mathit{\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 ${\mathfrak{L}}_{\mathrm{d}{x}_t}q$ is calculated using Cartan's formula:

$$\begin{array}{}\text{(64)}& {\mathfrak{L}}_{\mathrm{d}{x}_t}q=\mathrm{d}\left(i_{\mathrm{d}{x}_t}q\right)+i_{\mathrm{d}{x}_t}\mathrm{d}q.\end{array}$$

Essentially, the Lie derivative ${\mathfrak{L}}_{\mathrm{d}{x}_t}q$ corresponds to $T_t^{*}q-q+f^q\left(S\right)\mathrm{d}t$, if we assume that the deterministic forecast of q is simply the advection of q by u. More generally, ${\mathfrak{L}}_{\mathrm{d}{x}_t-u\mathrm{d}t}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 $\mathrm{m}=u^j\mathrm{d}{x}^j=u^{\mathrm{1}}\mathrm{d}{x}^{\mathrm{1}}+\mathrm{\dots }+u^n\mathrm{d}{x}^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 ${\mathfrak{L}}_{\mathrm{d}{x}_t}\mathit{\theta }$, where θ is a differential 1 form related to the momentum 1 form. For instance, θ=m in the example of “Stratonovich stochastic Euler–Poincaré flow” in Holm (2015), and $\mathit{\theta }=\mathrm{m}+R^j\mathrm{d}{x}^j$ in the example of “stochastic Euler–Boussinesq equations of a rotating stratified incompressible fluid” in Holm (2015). Already pointed out, the operator ${\mathfrak{L}}_{\mathrm{d}{x}_t}$ is closely related to $T_t^{*}$, and the momentum equation in SALT can be derived using our proposed scheme by properly choosing θ.

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 ${\mathfrak{L}}_{\mathrm{d}{x}_t}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 ${\mathfrak{L}}_{\mathrm{d}{x}_t-u\mathrm{d}t}$ 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)

$$\begin{array}{}\text{(65)}& {\displaystyle }{\partial }_{t}f+{\mathit{w}}^{\star }\cdot \mathrm{\nabla }f=\mathrm{\nabla }\cdot \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{a}\mathrm{\nabla }f\right)-\mathit{\sigma }\dot{\mathit{B}}\cdot \mathrm{\nabla }f\text{(66)}& {\displaystyle }{\mathit{w}}^{\star }=\mathit{w}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\nabla }\cdot \mathit{a}{)}^{\top }+\mathit{\sigma }(\mathrm{\nabla }\cdot \mathit{\sigma }{)}^{\top },\end{array}$$

where f can be any quantity that is assumed to be transported by the flow, i.e., $Df/Dt=\mathrm{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_idη_i is denoted $\mathit{\sigma }\mathrm{d}\mathit{B}={\mathit{\sigma }}_{•i}\mathrm{d}{B}_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

$$\begin{array}{}\text{(67)}& {\displaystyle }{\partial }_{t}f+{\mathit{w}}_S\cdot \mathrm{\nabla }f={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{•i}\cdot \mathrm{\nabla })({\mathit{\sigma }}_{•i}\cdot \mathrm{\nabla }f)-\left(\mathit{\sigma }\dot{\mathit{B}}\right)\cdot \mathrm{\nabla }f,\text{(68)}& {\displaystyle }=-(\mathit{\sigma }\circ \dot{\mathit{B}})\cdot \mathrm{\nabla }f,\text{(69)}& {\displaystyle }{\mathit{w}}_S=\mathit{w}+{\mathit{w}}_S^c\text{(70)}& {\displaystyle }{\mathit{w}}_S^c=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\nabla }\cdot \mathit{a}{)}^{\top }+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\sigma }(\mathrm{\nabla }\cdot \mathit{\sigma }{)}^{\top },\text{(71)}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{•i}\cdot \mathrm{\nabla }){\mathit{\sigma }}_{•i},\end{array}$$

where $\mathit{\sigma }\circ \dot{\mathit{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

$$\begin{array}{}\text{(72)}& {\mathrm{d}}^{\text{LU}}f=g^f\left(S\right)\mathrm{d}t+{\mathrm{d}}_s^{\text{LU}}f,\end{array}$$

where

$$\begin{array}{}\text{(73)}& {\displaystyle }{g}^f\left(S\right)=-\mathit{w}\cdot \mathrm{\nabla }f\text{(74)}& {\displaystyle }\begin{array}{rl}{\mathrm{d}}_s^{\text{LU}}f& =-{\mathit{w}}_S^c\cdot \mathrm{\nabla }f\mathrm{d}t+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{•i}\cdot \mathrm{\nabla })({\mathit{\sigma }}_{•i}\cdot \mathrm{\nabla }f)\mathrm{d}t\\ & -\left(\mathit{\sigma }\mathrm{d}\mathit{B}\right)\cdot \mathrm{\nabla }f.\end{array}\end{array}$$

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

$$\begin{array}{rl}-{\mathit{w}}_S^c\cdot \mathrm{\nabla }f\mathrm{d}t=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^q{\partial }_{x_q}{e}_i^p{\partial }_{x_p}f\\ {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\mathit{\sigma }}_{•i}\cdot \mathrm{\nabla })({\mathit{\sigma }}_{•i}\cdot \mathrm{\nabla }f)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{\partial }_{x_p}\left(e_i^q{\partial }_{x_q}f\right)\\ =& {\displaystyle \frac{\mathrm{1}}{\mathrm{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)\\ -\mathit{\sigma }\mathrm{d}\mathit{B}\cdot \mathrm{\nabla }f=& e_i^p{\partial }_{x_p}f\mathrm{d}{\mathit{\eta }}_i.\end{array}$$

Hence,

$$\begin{array}{}\text{(75)}& \begin{array}{rl}{\mathrm{d}}_s^{\text{LU}}f& =(e_i^q{\partial }_{x_q}{e}_i^p{\partial }_{x_p}f+{\displaystyle \frac{\mathrm{1}}{\mathrm{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}{\mathit{\eta }}_i.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(76)}& \begin{array}{rl}{T}_t\left(x\right)& =x+e_i^q{\partial }_{x_q}{e}_i\mathrm{\Delta }t+e_i\mathrm{\Delta }{\mathit{\eta }}_i=x-{\mathit{w}}_S^c\mathrm{\Delta }t\\ & +(-{\mathit{w}}_S^c\mathrm{\Delta }t-\mathit{\sigma }\mathrm{\Delta }\mathit{B}).\end{array}\end{array}$$

The LU equation can thus be derived by choosing θ=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 $(-{\mathit{w}}_S^c\mathrm{\Delta }t-\mathit{\sigma }\mathrm{\Delta }\mathit{B})=(\frac{\mathrm{1}}{\mathrm{2}}{e}_i^q{\partial }_{x_q}{e}_i\mathrm{\Delta }t+e_i\mathrm{\Delta }{\mathit{\eta }}_i)$ is the Itō noise and its Itō-to-Stratonovich correction. Hence, it corresponds to the Stratonovich noise e_i∘dη_i of the flow associated with T_t. The additional drift $-{\mathit{w}}_S^c\mathrm{\Delta }t$ is different in nature. It is related to the advection correction ${\mathit{w}}_S^c\cdot \mathrm{\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 ${\mathit{w}}_S=\mathit{w}+{\mathit{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),

$$\begin{array}{}\text{(77)}& T_t^{-\mathrm{1}}\left(x\right)=x-e_i\mathrm{\Delta }{\mathit{\eta }}_i=x+\mathit{\sigma }\mathrm{\Delta }\mathit{B}.\end{array}$$

Considering T_t to represent how much the model forecast differs from the true forecast at every time step, $T_t^{-\mathrm{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 θ=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 $\mathit{\rho }\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$ and $\mathit{\rho }\mathit{w}\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^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 $\mathit{\theta }=f\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^n$. Assuming an integral conservation $\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{V\left(t\right)}f=\mathrm{0}$ on a spatial domain V(t) transported by the flow, that theorem leads to the following SPDE:

$$\begin{array}{}\text{(78)}& \begin{array}{rl}{\displaystyle \frac{Df}{Dt}}& +\mathrm{\nabla }\cdot ({\mathit{w}}^{\star }+\mathit{\sigma }\dot{\mathit{B}})f={\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}〈\underset{\mathrm{0}}{\overset{t}{\int }}{D}_{t}f,\underset{\mathrm{0}}{\overset{t}{\int }}\mathrm{\nabla }\cdot \mathit{\sigma }\dot{\mathit{B}}〉\\ & =(\mathrm{\nabla }\cdot {\mathit{\sigma }}_{•i})(\mathrm{\nabla }\cdot {\mathit{\sigma }}_{•i}{)}^{T}f,\end{array}\end{array}$$

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):

$$\begin{array}{}\text{(79)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}f& +\mathrm{\nabla }\cdot \left({\mathit{w}}_{S}f\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\nabla }\cdot \left(\mathit{a}\mathrm{\nabla }f\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\nabla }\cdot \left({\mathit{\sigma }}_{•i}\right(\mathrm{\nabla }\cdot {\mathit{\sigma }}_{•i}{)}^{T}f)\\ & -\mathrm{\nabla }\cdot \left(\mathit{\sigma }\dot{\mathit{B}}f\right)\end{array}\text{(80)}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\nabla }\cdot \left({\mathit{\sigma }}_{•i}\right(\mathrm{\nabla }\cdot \left({\mathit{\sigma }}_{•i}f\right){)}^T)-\mathrm{\nabla }\cdot (\mathit{\sigma }\dot{\mathit{B}}f)\text{(81)}& {\displaystyle }=-\mathrm{\nabla }\cdot (\mathit{\sigma }\circ \dot{\mathit{B}}f).\end{array}$$

The original deterministic equation and stochastic perturbation correspond to

$$\begin{array}{}\text{(82)}& {\displaystyle }{g}^f\left(S\right)=-\mathrm{\nabla }\cdot \left(\mathit{w}f\right)\text{(83)}& {\displaystyle }\begin{array}{rl}{\mathrm{d}}_s^{\text{LU}}f& =(-\mathrm{\nabla }\cdot ({\mathit{w}}_S^{c}f)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\nabla }\cdot (\mathit{a}\mathrm{\nabla }f)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\nabla }\\ & \cdot \left({\mathit{\sigma }}_{•i}\right(\mathrm{\nabla }\cdot {\mathit{\sigma }}_{•i}{)}^{T}f\left)\right)\mathrm{d}t-\mathrm{\nabla }\cdot \left(\mathit{\sigma }\mathrm{d}\mathit{B}f\right)\end{array}\text{(84)}& {\displaystyle }=\mathrm{\nabla }\cdot \left(\left({\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\nabla }\cdot \mathit{a}\right)}^T\mathrm{d}t-\mathit{\sigma }\mathrm{d}\mathit{B}\right)f\right)+\mathrm{\nabla }\cdot \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{a}\mathrm{\nabla }f\right)\mathrm{d}t.\end{array}$$

Identifying $\mathit{a}={\mathit{\sigma }}_{•i}{\mathit{\sigma }}_{•i}^T=e_i{e}_i^T$ and $\mathit{\sigma }\dot{\mathit{B}}=-e_{i}d{\mathit{\eta }}_i$, Eq. (35) corresponds to Example 3.2.2 about n forms, with

$$\begin{array}{}\text{(85)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{V}& =-a^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\partial }_{x^q}{e}_i^p{e}_i^q-e_i^p{\partial }_{x^q}{e}_i^q)-e_i^p{\displaystyle \frac{\mathrm{d}{\mathit{\eta }}_i}{\mathrm{d}t}}\\ & =-{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{\nabla }\cdot \mathit{a}\right)}^T+\mathit{\sigma }\dot{\mathit{B}};\end{array}\end{array}$$

i.e.,

$$\begin{array}{}\text{(86)}& a^p={\partial }_{x^q}\left(e_i^p{e}_i^q\right)-\left(e_i^p{\partial }_{x^q}{e}_i^q\right)=e_i^q{\partial }_{x^q}{e}_i^p.\end{array}$$

Again, the remapping is obtained,

$$\begin{array}{}\text{(87)}& \begin{array}{rl}{T}_t\left(x\right)& =x+e_i^q{\partial }_{x_q}{e}_i\mathrm{\Delta }t+e_i\mathrm{\Delta }{\mathit{\eta }}_i=x-{\mathit{w}}_S^c\mathrm{\Delta }t\\ & +(-{\mathit{w}}_S^c\mathrm{\Delta }t-\mathit{\sigma }\mathrm{\Delta }\mathit{B}),\end{array}\end{array}$$

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{array}{}\text{(88)}& \left\{\begin{array}{l}\mathrm{\nabla }\cdot \mathit{\sigma }=\mathrm{0}\\ \mathrm{\nabla }\cdot \mathrm{\nabla }\cdot \mathit{a}=\mathrm{0}.\end{array}\right.\end{array}$$

Translating it into our notation, it reads

$$\left\{\begin{array}{l}{\partial }_{x_p}{e}_i^p=\mathrm{0}\text{ for each}i\\ {\partial }_{x_p}{\partial }_{x_q}\left(e_i^p{e}_i^q\right)=\mathrm{0}.\end{array}\right.$$

Applying the result in Example 3.1.2, straightforward calculation means Eq. (88) is equivalent to $T_t^{*}\mathit{\theta }=\mathit{\theta }$ for $\mathit{\theta }=\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^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):

$$\begin{array}{}\text{(89)}& {\displaystyle \frac{\partial h}{\partial t}}+\mathrm{\nabla }\cdot \left(h\stackrel{\mathrm{‾}}{u}\right)=\mathrm{0},\text{(90)}& {\displaystyle \frac{\partial \mathrm{\Theta }}{\partial t}}+(\stackrel{\mathrm{‾}}{u}\cdot \mathrm{\nabla })\mathrm{\Theta }=-\mathit{\kappa }(h\mathrm{\Theta }-h_{\mathrm{0}}{\mathrm{\Theta }}_{\mathrm{0}}),\text{(91)}& {\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u}}{\partial t}}+(\stackrel{\mathrm{‾}}{u}\cdot \mathrm{\nabla })\stackrel{\mathrm{‾}}{u}+f\widehat{z}\times \stackrel{\mathrm{‾}}{u}=-\mathrm{\nabla }\left(h\mathrm{\Theta }\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}h\mathrm{\nabla }\mathrm{\Theta }.\end{array}$$

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 ρ(x,t), while the lower layer is quiescent with a fixed constant density ρ₀. The state variable h represents the height of the active layer, and $\mathrm{\Theta }=g({\mathit{\rho }}_{\mathrm{0}}-\mathit{\rho })/{\mathit{\rho }}_{\mathrm{0}}$ is the density contrast. $\stackrel{\mathrm{‾}}{u}$ is the averaged horizontal velocity of the active layer at each column. Note that ρ<ρ₀ (hence Θ>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.

$$\begin{array}{}\text{(92)}& {\displaystyle }\text{Total energy: }E=\underset{\mathrm{\Omega }}{\int }{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(h\right|\stackrel{\mathrm{‾}}{u}{|}^{\mathrm{2}}+h^{\mathrm{2}}\mathrm{\Theta })d^{\mathrm{2}}x.\text{(93)}& {\displaystyle }\text{Total mass: }\mathcal{M}=\underset{\mathrm{\Omega }}{\int }h{d}^{\mathrm{2}}x.\text{(94)}& {\displaystyle }\text{Total momentum: }\mathrm{M}=\underset{\mathrm{\Omega }}{\int }h\stackrel{\mathrm{‾}}{u}{d}^{\mathrm{2}}x.\end{array}$$

The objective is thus to choose proper tensor fields ${\mathit{\theta }}_{\stackrel{\mathrm{‾}}{u}},{\mathit{\theta }}_h$, and θ_Θ for the state variables $\stackrel{\mathrm{‾}}{u},h,$ and Θ, 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 θ_h is ${\mathit{\theta }}_h=h\mathrm{d}{x}^{\mathrm{1}}\wedge \mathrm{d}{x}^{\mathrm{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 ${\mathit{\theta }}_{\stackrel{\mathrm{‾}}{u}}$ to be a function (differential 0 form). Therefore, the only choice for ${\mathit{\theta }}_{\stackrel{\mathrm{‾}}{u}}$ is ${\mathit{\theta }}_{\stackrel{\mathrm{‾}}{u}}=\stackrel{\mathrm{‾}}{u}$. This choice of ${\mathit{\theta }}_{\stackrel{\mathrm{‾}}{u}}$ and θ_h implies that $h|\stackrel{\mathrm{‾}}{u}{|}^{\mathrm{2}}$ also corresponds to a 2 form $|\stackrel{\mathrm{‾}}{u}{|}^{\mathrm{2}}{\mathit{\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 θ_Θ so that h²Θ corresponds to a differential 2 form. Note that θ_h is already a 2 form. We must thus select θ_Θ so that hΘ corresponds to a function. The only choice for θ_Θ is the contravariant tensor ${\mathit{\theta }}_{\mathrm{\Theta }}=\mathrm{\Theta }\frac{\partial }{\partial x^{\mathrm{1}}}\wedge \frac{\partial }{\partial x^{\mathrm{2}}}$. In this case, hΘ corresponds to the differential 0 form 〈θ_hθ_Θ〉=hΘ, where 〈,〉 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:

$$\begin{array}{}\text{(95)}& {\displaystyle }{\mathit{\theta }}_h=hd{x}^{\mathrm{1}}\wedge \mathrm{d}{x}^{\mathrm{2}}\text{(96)}& {\displaystyle }{\mathit{\theta }}_{{\stackrel{\mathrm{‾}}{u}}^j}={\stackrel{\mathrm{‾}}{u}}^j(\text{for}j=\mathrm{1},\mathrm{2})\text{(97)}& {\displaystyle }{\mathit{\theta }}_{\mathrm{\Theta }}=\mathrm{\Theta }{\displaystyle \frac{\partial }{\partial x^{\mathrm{1}}}}\wedge {\displaystyle \frac{\partial }{\partial x^{\mathrm{2}}}}.\end{array}$$

For

$$\begin{array}{}\text{(98)}& T_t\left(x\right)=x+a\mathrm{\Delta }t+e_i\mathrm{\Delta }{\mathit{\eta }}_i,\end{array}$$

we have

$$\begin{array}{}\text{(99)}& T_t^{-\mathrm{1}}\left(x\right)=x+(-a+e_i^p{\partial }_{x_p}{e}_i)\mathrm{\Delta }t-e_i\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}$$

Then $T_t^{*}{\mathit{\theta }}_h$, $T_t^{*}{\mathit{\theta }}_{\stackrel{\mathrm{‾}}{u}}$, and $(T_t^{-\mathrm{1}}{)}_{*}{\mathit{\theta }}_{\mathrm{\Theta }}$ can be calculated following Examples 3.1.3, 3.1.1, and 3.1.5. This further implies ${\mathrm{d}}_{s}h,{\mathrm{d}}_s\stackrel{\mathrm{‾}}{u}$, and d_sΘ, as shown in Examples 3.2.2, 3.2.1, and 3.2.4. Note that $T_t^{-\mathrm{1}}$ instead of T_t is applied to θ_Θ, as shown in Eq. (7). Finally, we end up with the following SPDE:

$$\begin{array}{}\text{(100)}& {\displaystyle }\begin{array}{rl}\mathrm{d}h=& -\mathrm{\nabla }\left(h\stackrel{\mathrm{‾}}{u}\right)\mathrm{d}t+\left(h\right({\partial }_{x_p}{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_i)+a^p{\partial }_{x_p}h\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{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)\mathrm{d}t\\ & +(h{\partial }_{x_p}{e}_i^p+{\partial }_{x_p}h{e}_i^p)\mathrm{d}{\mathit{\eta }}_i\end{array}\text{(101)}& {\displaystyle }\begin{array}{rl}\mathrm{d}\mathrm{\Theta }& =\mathit{\{}-(\stackrel{\mathrm{‾}}{u}\cdot \mathrm{\nabla })\mathrm{\Theta }-\mathit{\kappa }(h\mathrm{\Theta }-h_{\mathrm{0}}{\mathrm{\Theta }}_{\mathrm{0}}\left)\mathit{\right\}}\mathrm{d}t\\ & +\left(\mathrm{\Theta }\right(-{\partial }_{x_p}{a}^p+{\partial }_{x_p}({\partial }_{x_q}{e}_i{e}_i^q{)}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{J}_i)\\ & +{\partial }_{x_p}\mathrm{\Theta }{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x_p}{\partial }_{x_q}\mathrm{\Theta }-{\partial }_{x_p}\mathrm{\Theta }{e}_i^p{\partial }_{x_q}{e}_i^q)\mathrm{d}t\\ & -(\mathrm{\Theta }{\partial }_{x_p}{e}_i^p-{\partial }_{x_p}\mathrm{\Theta }{e}_i^p)\mathrm{d}{\mathit{\eta }}_i\end{array}\text{(102)}& {\displaystyle }\begin{array}{rl}\mathrm{d}{\stackrel{\mathrm{‾}}{u}}^j=& -\mathit{\left\{}\right(\stackrel{\mathrm{‾}}{u}\cdot \mathrm{\nabla })\stackrel{\mathrm{‾}}{u}-f\widehat{z}\times \stackrel{\mathrm{‾}}{u}-\mathrm{\nabla }(h\mathrm{\Theta })+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}h\mathrm{\nabla }\mathrm{\Theta }{\mathit{\}}}^j\mathrm{d}t\\ & +({\partial }_{x_p}{\stackrel{\mathrm{‾}}{u}}^j{a}^p+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^p{e}_i^q{\partial }_{x_p}{\partial }_{x_q}{\stackrel{\mathrm{‾}}{u}}^j)\mathrm{d}t\\ & +{\partial }_{x_p}{\stackrel{\mathrm{‾}}{u}}^j{e}_i^p\mathrm{d}{\mathit{\eta }}_i,\end{array}\end{array}$$

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 θ by $\mathit{\theta }\leftarrow T_t^{*}\mathit{\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

$$\begin{array}{}\text{(103)}& \mathrm{d}S=f\left(S\right)\mathrm{d}t+{\mathrm{d}}_{s}S,\end{array}$$

where f(S)dt is the incremental of S given by the original deterministic system. The term d_sS 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 θ. A key point is indeed to link the state variable S with some tensor field θ. The choice of θ 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 θ. 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 θ 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 ${\mathfrak{L}}_{\mathrm{d}{x}_t-u\mathrm{d}t}$, 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^{-\mathrm{1}}$

Suppose that

$$\begin{array}{}\text{(A1)}& T_t\left(x\right)=x+a\mathrm{\Delta }t+e_i\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}$$

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

$$\begin{array}{}\text{(A2)}& T_t^{-\mathrm{1}}\left(x\right)=x+z\mathrm{\Delta }t+b_i\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}$$

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

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}x=& T_t\left(T_t^{-\mathrm{1}}\right(x\left)\right)=T_t(x+z\mathrm{\Delta }t+b_i\mathrm{\Delta }{\mathit{\eta }}_i)\\ =& x+z\mathrm{\Delta }t+b_i\mathrm{\Delta }{\mathit{\eta }}_i+a{|}_{x+z\mathrm{\Delta }t+b_i\mathrm{\Delta }{\mathit{\eta }}_i}\mathrm{\Delta }t\\ & +e_i{|}_{x+z\mathrm{\Delta }t+b_i\mathrm{\Delta }{\mathit{\eta }}_i}\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}\end{array}$$

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{array}{}\text{(A4)}& \begin{array}{rl}a{|}_{x+a\mathrm{\Delta }t+b_i\mathrm{\Delta }{\mathit{\eta }}_i}\mathrm{\Delta }t=& a{|}_x\mathrm{\Delta }t+o\left(\mathrm{\Delta }t\right)\\ e_i{|}_{x+z\mathrm{\Delta }t+b_i\mathrm{\Delta }{\mathit{\eta }}_i}\mathrm{\Delta }{\mathit{\eta }}_i=& e_i{|}_x\mathrm{\Delta }{\mathit{\eta }}_i+e_{ip}{b}_i^p{|}_x\mathrm{\Delta }t+o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

Therefore,

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}x& =T_t\left(T_t^{-\mathrm{1}}\right(x\left)\right)=x+(z+a+e_{ip}{b}_i^p)\mathrm{\Delta }t\\ & +(b_i+e_i)\mathrm{\Delta }{\mathit{\eta }}_i+o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

This implies that

$$\begin{array}{}\text{(A6)}& {\displaystyle }{b}_i+e_i=\mathrm{0}\text{(A7)}& {\displaystyle }z+a+e_{ip}{b}_i^p=\mathrm{0}.\end{array}$$

Therefore,

$$\begin{array}{}\text{(A8)}& {\displaystyle }{b}_i=-e_i\text{(A9)}& {\displaystyle }z=-a+e_{ip}{e}_i^p,\end{array}$$

or equivalently,

$$\begin{array}{}\text{(A10)}& T_t^{-\mathrm{1}}\left(x\right)=x+(-a+e_{ip}{e}_i^p)\mathrm{\Delta }t-e_i\mathrm{\Delta }{\mathit{\eta }}_i.\end{array}$$

Appendix B: Derivation of $T_t^{*}\mathit{\theta }$

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

$$\begin{array}{}\text{(B1)}& \mathit{\theta }=\sum _{i_{\mathrm{1}}<\mathrm{\dots }<i_k}{f}^{i_{\mathrm{1}},\mathrm{\dots },i_k}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}.\end{array}$$

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

$$\begin{array}{}\text{(B2)}& \mathit{\theta }=f\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}\end{array}$$

for some $\mathrm{1}\le i_{\mathrm{1}}<\mathrm{\dots }<i_k\le n$. Let $T_t\left(x\right)=\left(T_t^{\mathrm{1}}\right(x),\mathrm{\dots },T_t^n(x\left)\right)$, then

$$\begin{array}{}\text{(B3)}& \left(T_t^{*}\mathit{\theta }\right)\left(x\right)=f\left(T_t\right(x\left)\right)\mathrm{d}{T}_t^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{T}_t^{i_k}.\end{array}$$

We calculate f(T_t(x)) and $\mathrm{d}{T}_t^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{T}_t^{i_k}$ separately. We denote $\mathrm{\Delta }x=T_t\left(x\right)-x=a\mathrm{\Delta }t+e_i\mathrm{\Delta }{\mathit{\eta }}_i$ and H_f the Hessian matrix of f. At a given time t, f is assumed independent from the noise Δη_i(t). Then

$$\begin{array}{}\text{(B4)}& {\displaystyle }\begin{array}{rl}f\left(T_t\right(x\left)\right)& =f(x+\mathrm{\Delta }x)=f\left(x\right)+\langle \mathrm{\nabla }f\mathrm{\Delta }x\rangle +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{\Delta }x{)}^{\top }\\ & H_f\mathrm{\Delta }x+o\left(\right(\mathrm{\Delta }x{)}^{\mathrm{2}})\end{array}\text{(B5)}& {\displaystyle }=f\left(x\right)+\langle \mathrm{\nabla }fa\mathrm{\Delta }t+e_i\mathrm{\Delta }{\mathit{\eta }}_i\rangle +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^{\top }{H}_f{e}_i(\mathrm{\Delta }{\mathit{\eta }}_i{)}^{\mathrm{2}}\text{(B6)}& {\displaystyle }\begin{array}{rl}+\mathcal{O}& \left(\right(\mathrm{\Delta }t{)}^{\mathrm{2}})+\mathcal{O}(\mathrm{\Delta }t\mathrm{\Delta }{\mathit{\eta }}_i)+o(\left(\mathrm{\Delta }t{)}^{\mathrm{2}}\right)+o\left(\right(\mathrm{\Delta }{\mathit{\eta }}_i{)}^{\mathrm{2}})\\ & +o\left(\mathrm{\Delta }t\mathrm{\Delta }{\mathit{\eta }}_i\right).\end{array}\end{array}$$

According to Itô's lemma, dηdη=dt, and we can replace (Δη_i)² with Δt. Hence,

$$\begin{array}{}\text{(B7)}& {\displaystyle }\begin{array}{rl}f\left(T_t\right(x\left)\right)& =f\left(x\right)+\langle \mathrm{\nabla }fa\rangle \mathrm{\Delta }t+\langle \mathrm{\nabla }f{e}_i\rangle \mathrm{\Delta }{\mathit{\eta }}_i\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^{\top }{H}_f{e}_i\mathrm{\Delta }t+o\left(\mathrm{\Delta }t\right)\end{array}\text{(B8)}& {\displaystyle }\begin{array}{rl}=f\left(x\right)& +(\langle \mathrm{\nabla }fa\rangle +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i{H}_f{e}_i)\mathrm{\Delta }t+\langle \mathrm{\nabla }f{e}_i\rangle \mathrm{\Delta }{\mathit{\eta }}_i\\ & +o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

Next,

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}& T_t^{*}(\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})=\mathrm{d}{T}_t^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{T}_t^{i_k}\\ =& (\mathrm{d}{x}^{i_{\mathrm{1}}}+d{a}^{i_{\mathrm{1}}}\mathrm{\Delta }t+d{e}_i^{i_{\mathrm{1}}}\mathrm{\Delta }{\mathit{\eta }}_i)\wedge \mathrm{\dots }\wedge \\ & (\mathrm{d}{x}^{i_k}+d{a}^{i_k}\mathrm{\Delta }t+d{e}_i^{i_k}\mathrm{\Delta }{\mathit{\eta }}_i).\end{array}\end{array}$$

Note that $d{a}^{i_j}$ and $d{e}_i^{i_j}$ refer to the spatial differentiation. Again, we apply the “discrete version” of Itô's rule (Δη_i)²=Δt and collect all the terms of order 𝒪(Δt) and 𝒪(Δη_i):

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}{T}_t^{*}& (\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})=\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}\\ & +(\sum _{s=\mathrm{1}}^k\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge d{a}^{i_s}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\mathrm{\Delta }t\\ & +(\sum _{s=\mathrm{1}}^{k}d{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge d{e}_i^{i_s}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\mathrm{\Delta }{\mathit{\eta }}_i\\ & +(\sum _{s<r}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge d{e}_i^{i_s}\wedge \mathrm{\dots }\wedge d{e}_i^{i_r}\\ & \wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\mathrm{\Delta }t\\ & +o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

According to the chain rule, $\mathrm{d}{a}^{i_s}={\partial }_{x^j}{a}^{i_s}\mathrm{d}{x}^j$, $\mathrm{d}{e}_i^{i_s}={\partial }_{x^j}{e}_i^{i_s}\mathrm{d}{x}^j$. Note that ${\partial }_{x^j}{e}_i^{i_s}$ refers to the i_sth component of ${\partial }_{x^j}{e}_i$, where ${\partial }_{x^j}{e}_i=\frac{\partial e_i}{\partial x^j}$ and e_i(x)∈ℝⁿ is the ith basis vector field of T_t. Hence,

$$\begin{array}{}\text{(B11)}& \begin{array}{rl}{T}_t^{*}& (\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\\ =& \mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}+(\sum _{s=\mathrm{1}}^k{\partial }_{x^j}{a}^{i_s}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^j\\ & \wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\mathrm{\Delta }t\\ & +(\sum _{s=\mathrm{1}}^k{\partial }_{x^j}{e}_i^{i_s}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^j\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\mathrm{\Delta }{\mathit{\eta }}_i\\ & +(\sum _{s<r}{\partial }_{x^j}{e}_i^{i_s}{\partial }_{x^l}{e}_i^{i_r}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^j\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^l\\ & \wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\mathrm{\Delta }t\\ & +o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

Combining Eqs. (B8) and (B11), with application of Itô's lemma, all terms of order o(Δt) are then removed, to obtain

$$\begin{array}{}\text{(B12)}& \begin{array}{rl}{T}_t^{*}\mathit{\theta }=& f\left(T_t\right(x\left)\right)T_t^{*}(\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\\ =& \mathit{\theta }+\mathit{\left\{}\right(\langle \mathrm{\nabla }fa\rangle +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{e}_i^{\top }{H}_f{e}_i)\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_n}\\ & +\sum _{s=\mathrm{1}}^{k}f{\partial }_{x^j}{a}^{i_s}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\mathrm{d}{x}^j\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}\\ & +(\sum _{s<r}f{\partial }_{x^j}{e}_i^{i_s}{\partial }_{x^l}{e}_i^{i_r}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^j\wedge \mathrm{\dots }\wedge \\ & \mathrm{d}{x}^l\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k})\\ & +(\sum _{s=\mathrm{1}}^k\langle \mathrm{\nabla }f{e}_i\rangle {\partial }_{x^j}{e}_i^{i_s}\mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^j\\ & \wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}\left)\mathit{\right\}}\mathrm{\Delta }t\\ & +\mathit{\{}\langle \mathrm{\nabla }f{e}_i\rangle \mathrm{d}{x}^{i_{\mathrm{1}}}\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}\\ & +\sum _{s=\mathrm{1}}^{k}f{\partial }_{x^j}{e}_i^{i_s}\mathrm{d}{x}^{i_{\mathrm{1}}}\\ & \wedge \mathrm{\dots }\wedge \mathrm{d}{x}^j\wedge \mathrm{\dots }\wedge \mathrm{d}{x}^{i_k}\mathit{\}}\mathrm{\Delta }{\mathit{\eta }}_i\\ & +o\left(\mathrm{\Delta }t\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B13)}& T_t^{*}\mathit{\theta }=\mathit{\theta }+\mathcal{M}\left(\mathit{\theta }\right)\mathrm{\Delta }t+{\mathcal{N}}_i\left(\mathit{\theta }\right)\mathrm{\Delta }{\mathit{\eta }}_i,\end{array}$$

for some differential k forms ℳ(θ) and 𝒩_i(θ).