Fluid flow dynamics under location uncertainty

# Fluid flow dynamics under location uncertainty

E. MÉMIN
INRIA
Corresponding author. Email: Etienne.Memin@inria.fr Inria
Fluminance group    Campus universitaire de Beaulieu    35042 Rennes Cedex    France
Received 11 February 2013; in final form 26 July 2013; first published online ????
###### Abstract

We present a derivation of a stochastic model of Navier Stokes equations that relies on a decomposition of the velocity fields into a differentiable drift component and a time uncorrelated uncertainty random term. This type of decomposition is reminiscent in spirit to the classical Reynolds decomposition. However, the random velocity fluctuations considered here are not differentiable with respect to time, and they must be handled through stochastic calculus. The dynamics associated with the differentiable drift component is derived from a stochastic version of the Reynolds transport theorem. It includes in its general form an uncertainty dependent subgrid bulk formula that cannot be immediately related to the usual Boussinesq eddy viscosity assumption constructed from thermal molecular agitation analogy. This formulation, emerging from uncertainties on the fluid parcels location, explains with another viewpoint some subgrid eddy diffusion models currently used in computational fluid dynamics or in geophysical sciences and paves the way for new large-scales flow modeling. We finally describe an applications of our formalism to the derivation of stochastic versions of the Shallow water equations or to the definition of reduced order dynamical systems.

Keywords: Damping; Porous structure; Reflection; Transmission; Matching conditions Stochastic Reynolds Transport Theorem; Stochastic Navier-Stokes equation; Large scales flow modeling; Uncertainty specification, Subgrid stress modeling

\jvol

00 \jnum00 \jyear2012

## 1 Introduction

For several years there has been a deep and growing interest in defining stochastic models for climate or geophysical sciences (Majda et al., 1999; Slingo and Palmer, 2011) – see also on this subject the thematic issue (Palmer and Williams, 2008) and references therein or the review (Frederiksen et al., 2013). This interest is strongly motivated by the necessity of devising large-scale evolution models that allow taking into account statistical descriptions of processes that cannot be accurately specified - to keep, for instance, an affordable computational time. This includes small-scale physical forcing (cloud convection, boundary layer turbulence, radiative interaction, etc.) and numerical hypothesis and choices operated at the dynamics modeling level, such as a scale coarsening in a given direction, or the selection of a particular functional space to express the solution. The modeling simplifications and the different unresolved processes involved introduce de facto errors and uncertainties on the constitutive equations of the state variables, which are, in general, so complex that only a probabilistic modeling can be envisaged.

In traditional large-scale modeling, the interaction between the resolved flow component and the unresolved components lies essentially in the constitution of a so-called subgrid stress tensor, which is usually not related to uncertainty or to an error concept. We believe that it is important to extend this notion in order to take into consideration in a more appropriate way the action, on the resolved component, of stochastic processes modeling errors and uncertainties of the state variables’ dynamics.

The modeling and the handling along time of such uncertainties are crucial, for instance, for ensemble forecasting issues in meteorology and oceanography, where ensembles of runs are generated through randomization of the dynamics’ parameters, accompanied eventually by stochastic forcing that mimics the effect of unresolved physical processes (Slingo and Palmer, 2011). Problems related to the underestimation of the ensemble spread and a lack of representativity of the subgrid stress in terms of the model errors constitute in particular problematic limitations of the ensemble techniques’ forecasting skill. Such stochastic evolution models are also needed for data assimilation procedures defined through Monte Carlo implementation of stochastic filters referred to as ensemble filters or particle filters in the literature (Evensen, 2006; Gordon et al., 2001). In all these situations, a large-scale flow evolution description that includes stochastic forcing terms and an uncertainty related subgrid stress expression hence appears necessary. The root problem of this general issue lies essentially in the construction of adequate stochastic versions of the Navier-Stokes equations.

The establishment of sound stochastic dynamical models to describe the fluid flow but also the evolution of the different random terms encoding uncertainties or errors is a difficult issue compounded with an involved mathematical analysis. As initiated by (Bensoussan and Temam, 1973) and intensively studied by many authors since then (see for instance the review Flandoli, 2008, and references therein) such a formalization has been mainly considered through the addition of random forcing terms to a standard expression of Navier-Stokes equations. However, this construction is limited by a priori assumptions about the noise structure. Furthermore, the question whether this noise should be multiplicative or additive immediately arises.

Another family of methods initiated by Kraichnan (1959) consists to close the large scale representation by neglecting statistical correlations in the Fourier space through the so-called Direct-interaction approximation (DIA). These methods can be generalized by relying on a Langevin stochastic representation (Kraichnan, 1970; Leith, 1971) to devise advanced stochastic subgrid models for barotropic flows or quasi-geostrophic models where interactions between turbulent eddies and topography are taken into account(Frederiksen, 1999, 2012; Frederiksen et al., 2013). This strategy based on renormalized perturbation theory comes to randomize the spectral Navier-Sokes representation by replacing the nonlinear interaction terms by appropriate random forcing and damping terms.

Compared to these works, we wish to explore the problem via a somewhat reverse strategy. Instead of considering a given – eventually simplified – dynamics and then to supplement it with random forcing terms to model errors carried by unresolved or unknown processes, we will start from a general Lagrangian stochastic description that incorporates uncertainties on the fluid motion. The sought-after Eulerian dynamics is then deduced from this general stochastic velocity description and standard physical principles or approximations. This construction, reminiscent to the framework proposed by (Mikulevicius and Rozovskii, 2004) and initiated by (Brzeźniak et al., 1991), has the great advantage to let naturally emerge deterministic and stochastic uncertainty terms related to the different errors transported by the evolution model. Deterministic approximations or stochastic implementations of the dynamics can then be considered on solid grounds.

Following this route, we aim in this paper at devising stochastic dynamics for the description of fluid flows under uncertainties. Such uncertainties, modeled through the introduction of random terms, allow taking into account approximations or truncation effects performed during the dynamics constitution steps. This includes for instance the modeling of the unresolved scales interaction in Large-Eddies Simulations (LES) or Reynolds Average Numerical Simulations (RANS). These uncertainties may also encode numerical errors, unknown forcing terms or uncertainties on the initial conditions. They gather the effects of processes one does not wish to accurately model. However, as they propagate along time and interact with the resolved components their effects must be properly taken into account.

We will assume throughout this study that the fluid particles displacements can be separated in two components: a smooth differentiable function of time and space and an uncertainty function uncorrelated in time but correlated in space. This latter component is formulated as a function of Brownian motion and the whole displacement is defined as an Ito diffusion of the form:

 dX(x,t)=w(X(x,t),t)dt+σ(X(x,t),t)d^Bt, (1)

where is the fluid flow map, which represents the trajectory followed by a fluid particle starting at point of the bounded domain . This constitutes a Lagrangian representation of the fluid dynamics and figures the Lagrangian displacement map of the flow at time . In expression (1), the velocity vector field, , corresponds to the smooth resolved velocity component of the flow. It is assumed to be a deterministic function (of random arguments) and to have twice differentiable components: . When it does not depend on random arguments this velocity field represents the expectation of the whole random velocity field. This is the situation encountered in the case of mean field dynamics. The second term is a generalized random field that assembles the unresolved flow component and all the uncertainties we have on the flow. The combination of both velocity fields provides an Eulerian description of the complete velocity fields driving the particles:

 U(x,t)=w(x,t)dt+σ(x,t)d^Bt. (2)

The whole random field, , should be a solution of the Navier-Stokes equations and is defined here as the combination of a stochastic uncertainty component and a deterministic “resolved” component driven respectively by unknown characteristics, and that have to be determined or specified.

This framework is also related to the work of (Constantin and Iyer, 2008). Nevertheless, their study aimed at building a Monte-Carlo Lagrangian representation of the Navier-Stokes equations. It is limited to a constant diffusion tensor and exhibits additional difficulties to deal with boundary conditions (Constantin and Iyer, 2011). In this work, we do not look for a stochastic Lagrangian representation of the deterministic Navier-Stokes equation. We seek instead, a deterministic smooth representation of the drift velocity component corresponding to a solution of a stochastic Navier-Stokes equation. The goals are thus in some way opposite. Our objective can be interpreted in terms of Large Eddies Simulation (LES) or Reynolds Average Numerical Simulation (RANS), as we aim at separating a resolved flow component from an unresolved one, respectively defined as a differentiable component and a time uncorrelated uncertainty component. From this prospect, the approach proposed here is closer in spirit to RANS techniques than to the LES paradigm as no spatial filtering is applied (see (Lesieur and Métais, 1996; Meneveau and Katz, 2000; Sagaut, 2005), for extensive reviews on the subject). Furthermore, let us point out that the uncertainty component lives at all the hydrodynamical scales. Thus, there is no spatial scale separation principle here but rather a temporal decomposition in terms of a highly oscillating process with no time differentiation property and a smooth differentiable component. This approach is also close in spirit to the separation in term of a ”coherent” component plus noise operated through adaptive wavelet basis (Farge et al., 2001, 1999). However, contrary to this approach relying on a Galerkin projection with an adaptive scale thresholding, our decomposition makes appear a diffusion tensor assembling the action of the unresolved uncertainty component on the resolved component.

More precisely, the resolved component of the Navier-Stokes equation we consider includes in its general form an anisotropic diffusion term that emerges due to the presence of the random uncertainty term. By analogy with the Reynolds decomposition and LES, we refer to this tensor as the subgrid stress tensor. However, it is important to outline that its construction differs completely from the Reynold stress definition. As a matter of fact, in our case the unresolved fluctuating component is a non differentiable random process, and stochastic calculus differentiation rules have to be used. We will see that in the general case the resulting subgrid stress cannot be immediately related to the usual eddy viscosity assumption, formulated in the nineteen century by (Boussinesq, 1877) from thermal molecular agitation analogy – also commonly referred to as Boussinesq assumption in the litterature – and that is still intensively used in the Large Eddies Simulation paradigm since the work of (Smagorinsky, 1963) and (Lilly, 1966).

The simpler form of eddy diffusion will be recovered only when the uncertainty will be confined to homogeneous random fields such as the Kraichnan random field (Kraichnan, 1968). For particular forms of the uncertainty random component, this formulation will explain with another viewpoint some subgrid eddy diffusion models proposed in computational fluid dynamics or in geophysical sciences. Following this route, an appropriate definition of uncertainty models adapted to given situations allows opening the way for new large scales flow modeling.

In a very similar way as in the deterministic case, the representation we propose will be determined from a stochastic version of the Reynolds transport theorem relying on a specific model of the uncertainty random field.

## 2 Construction of the uncertainty random field

Before presenting in detail the derivation we consider to built stochastic fluid flow evolution models, we need to specify the random component that will encode the uncertainty associated to the fluid parcels location. This random field, which lives on the bounded domain , relies on a Brownian motion field, refer hereafter to as Brownian avatar (as a smooth idealization of a dense Brownian map) defined on . This Brownian field is defined in the following section. It is built from a finite dimensional discrete set of standard Brownian variable and tends in law to a continuous limit. This construction will allow us to handle the spatial derivatives of the uncertainty field with simple calculus.

### 2.1 Brownian motion field avatar

This random field, denoted , is built from a finite dimensional discrete set of standard Brownian variables as

 ^Bnt(x)=1√nn∑i=1Bt(xi)φν(x−xi), (3)

where is a set of independent -dimensional (with or 3) standard Brownian motions centered on the points of a discrete grid and is a Gaussian function of standard deviation . It is immediate to check that is a zero mean Gaussian process, with uncorrelated increments. Its spatial covariance tensor is defined as

 E[^Bnt(x)^BnTt(y)] =tnId∑iφν(x−xi)φν(y−xi), (4)

and, for an infinite number of grid points, tends to

 Q=limn→∞E[^Bnt(x)^BnTt(y)]=tφ√2ν(x−y)Id, (5)

as

 limn→∞1n∑iφν(x−xi)φν(y−xi) =(4πν2)−d/2exp(−14ν2∥x−y∥2) (6a) =φ√2ν(x−y). (6b)

It can be checked the operator

 Qf(x)=∫RdQ(x,y)f(y)dy, (7)

is symmetric positive definite. In addition, in order to define a well defined Gaussian process the trace of this covariance must be bounded for any orthonormal basis of the associated function space. It can be checked this covariance operator is not of finite trace in . However it is well defined in a bigger space (that includes in particular the constant functions) such as the space of mean square integrable functions associated to the measure . As a matter of fact, denoting the inner product associated to the space , we have for any orthonormal basis of this space

 trQ\lx@stackrel△=∞∑i=1(Qei,ei)α =td∞∑i=1∫Rd∫Rdφ√2ν(x−y)ei(x)ei(y)dμ(y)dμ(x) (8a) ≤td∞∑i=1(∫Rdφ2ν(0)ei(y)dμ(y))(∫Rdφ2ν(0)ei(x)dμ(x)) (8b) ≤td∞∑i=1(φ2ν(0),ei)2α (8c) ≤td∫Rdφ22ν(0)dμ(x)≤C,whenα>d−1. (8d)

As is a Gaussian process, it tends thus in law to a zero mean continuous process in with the same limiting covariance 111Note the covariance operator would be also well defined in a periodic domain . In this case .. This limiting process will be denoted in a formal way through a convolution product

 Bt⋆φν(x)=∫RdBt(x′)φν(x−x′)dx′, (9)

and will be identified to in the following. Furthermore, it can checked that the covariance of has a finite trace. Indeed, for any orthonormal basis of , denoting by the associated inner product and its induced norm, we have

 trE[^Bnt(x)^BnTt(x)] (10a) =∞∑k=1E(^Bnt,ek)2=E∥∥^Bnt∥∥22, (10b)

where the Brownian avatar norm is given by

 E∥∥^Bnt∥∥22 (11)

The energy of the Brownian avatar hence does not depend on the number of grid points used for its construction but only on the standard deviation of the Gaussian smoothing function. The trace of the limiting covariance is thus provided by (11). At finite time horizon, when the standard deviation tends to infinity, the energy brought by the uncertainty tends to zero and the SDE (1) behaves almost like a deterministic evolution law. At the opposite, in the zero limit of the smoothing function standard deviation, the variance – or the energy – becomes unbounded and the Ito stochastic integral can be no more defined. For non zero standard deviation the covariance operator is of finite trace and the Brownian avatar tends to a well defined -Wiener process (Prato and Zabczyk, 1992) in . Let us point out that such Gaussian correlation is probably the most widely used kernel in Geo-statistics and Machine Learning despite its strong unrealistic smoothing characteristics. In our case this kernel will be used in combination with a diffusion tensor.

### 2.2 White noise avatar

The analogue of the white noise on the bounded domain is similarly defined in a formal way from the generalized function and, , a linear bounded deterministic symmetric operator of the space of mean square integrable functions with null condition outside the domain interior . We will furthermore assume that this operator is Hilbert-Schmidt (for any orthonormal basis of of the operator has a bounded operator norm: ).

 σ(x,t)d^Bt=∫Ωσt(x,y)d^Bt(y)dy. (12)

This operator is referred in the following to as the diffusion tensor. The limiting covariance tensor of the uncertainty component , (also called the two points correlation tensor) reads

 Q =limn→∞1ndtδ(t−s)n∑i=0σ(x,∙,t)⋆φν(xi)σ(∙,y,t)⋆φν(xi) =dtδ(t−s)σφν(x,t)σφν(y,t), (13)

where is a filtered version of along its first or second component since it is symmetric. As the diffusion tensor is Hilbert-Schmidt, we observe immediately through Young’s inequality that the covariance is of finite trace. It constitutes thus a well defined covariance operator and the resulting process is a -Wiener process in .

Processes of central importance in stochastic calculus are the quadratic variation and co-variation processes. They correspond to the variance and covariance of the stochastic increments along time. The quadratic co-variation process denoted as , (respectively the quadratic variation for ) is defined as the limit in probability over a partition of with , and a partition spacing , noted as and such that when :

 ⟨X,Y⟩t=Plim|δt|n→0n−1∑i=0(X(ti+1)−X(ti))(Y(ti+1)−Y(ti))T.

For Brownian motion these covariations can be easily computed and are given by the following rules , , where is a deterministic function and a scalar Brownian motion. In our case, the quadratic variation of the uncertainty component reads

 ⟨∫t0(σ(x,t)d^Bt)i,∫t0(σ(x,t)d^Bt)j⟩ = limn→∞1n∫t0n∑ℓ=0∑kσik (x,∙,s)⋆φν(xℓ)σkj(∙,x,s)⋆φν(xℓ)ds =∫t0∑kσikφν(x,s)σkjφν(x,s)ds \lx@stackrel△=∫t0aij(x,s)dsa.s. (14)

Its time derivative corresponds to the diagonal of the spatial covariance tensor. We can note that for homogeneous random fields these expressions simplify. As a matter of fact, in this case, the random component has a covariance tensor that is defined as

 Q=dtδ(t−s)σ(∙,t)⋆σ(∙,t)⋆φ√2ν(x−y). (15)

It remains thus homogeneous and its quadratic variation tends to a spatially constant matrix for an infinite number of grid points:

 limn→∞1nn∑i=0[σ(∙,t)⋆φν(x−xi)]2dt= dt∫Ω[σ(∙,t)⋆φν(x)]2dx = Σ(t)dta.s. (16)

Besides, we remark that spatial derivatives of the white noise avatar can be defined and manipulated in a simple manner:

 \upartialxi(σ(x)d^Bt)= limn→∞1√nn∑j=0\upartialxiσφν(x,xj)dBt(xj) = \upartialxiσφν(x)d^Bt.

As a direct consequence, we note that a divergence free random component is thus necessarily associated to a divergence free diffusion tensor. The derivative operator of the identity diffusion tensor is also still well defined for a nonzero standard deviation and acts solely on the white noise avatar Gaussian smoothing function.

### 2.3 Kraichnan turbulent model

Toys models pioneered by (Kraichnan, 1968) and intensively explored for the theoretical analysis of passive scalar turbulence (Gawedzky and Kupiainen, 1995; Kraichnan, 1968; Majda and Kramer, 1999) can be easily specified with such a model. The Kraichnan random field model

 dξζt=P⋆ψγκ⋆fζ⋆dBt, (17)

is formally defined from the divergence free projector , a power function for and a band-pass cut-off function, , on the so-called inertial range that lies at high Reynolds number between the short dissipative scale and the large integral scale at which the forcing takes place. The constant has a dimension of , which follows from dimensional analysis (the energy transfer rate and ). The spectral correlation tensor of is defined as

 ˆQ(k)ij=|k|−ζ−d(δij−kikj|k|2)(ˆψγκ)2,

and the variance (or time derivative of the quadratic variation process) of this homogeneous random field reads

 d⟨ξζt(x),ξζt(x)⟩ij=dtδij∫Ω[P⋆ψγκ⋆fζ(x)]2dxa.s. (18)

which is constant and given for a band-pass filter, , defined as a box filter as

 dtCζ(2π)dd−1d2πd/2Γ(d/2)ζ−1(Lζ−ℓζD)δij. (19)

The square root of this quantity represents the mean absolute value of the velocity field components. It blows up when the integral scale tends to infinity and exhibits a strong dependency on the largest length scale. This injection scale dominates the turbulent energy and corresponds to the distance beyond which velocities are uncorrelated.

We will see that such random fields, which are both homogeneous in space and uncorrelated in time, lead to a simple form of the dynamics drift term involving an intuitive eddy viscosity term defined from the uncertainty random field variance. However more general random fields will let rise a formulation of the Navier-Stokes equations in which emerges an anisotropic diffusion term that cannot be immediately related to the usual eddy viscosity assumption first formulated in the nineteen century by (Boussinesq, 1877)222referred sometimes for that reason as the Boussinesq assumption. and that is still predominantly used in the Large Eddies Simulation paradigm.

## 3 Stochastic Reynolds transport theorem

In a similar way as in the deterministic case, our derivation relies essentially on a stochastic version of the Reynolds transport theorem. The rate of change of a scalar function within a material fluid volume transported by (1) is

 d∫V(t)q(x,t)dx=∫V(t){dtq+[∇⋅(qw)+12∥∇⋅σ∥2q−∑i,j12\upartial2\upartialxi\upartialxj(aijq)|∇⋅σ=0]dt+∇⋅(qσd^Bt)}dx. (20)

In this expression the forth right-hand term must be computed considering the diffusion tensor is divergence free, and the tensor denotes the uncertainty variance defined in (14) and denotes the function increment at a fixed point for an infinitesimal time step.

To give some insights on how this expression is obtained, let us consider a sufficiently spatially regular enough scalar function of compact support, transported by the flow and that vanishes outside and on . As this function is assumed to be transported by the stochastic flow (1), it constitutes a stochastic process defined from its initial time value :

 ϕ(Xt,t)=g(x0),

where both functions have bounded spatial gradients. The initial function vanishes outside the initial volume and on its boundary. Let us stress the fact that function cannot be a deterministic function. As a matter of fact, if it was the case, its differential would be given by a standard Ito formula

 dϕ(t,Xt)=(\upartialtϕ+∇ϕ⋅w+12∑i,jd⟨Xit,Xjt⟩\upartial2ϕ\upartialxi\upartialxj)dt+∇ϕ⋅σd^Bt.

This quantity cancels as is transported by the flow. A separation between the slow deterministic terms and the fast Brownian term would yield a very specific uncertainty lying on the function iso-values surfaces () or a null uncertainty, which is not the general goal followed here.

As is a random function, the differential of involves the composition of two stochastic processes. Its evaluation requires the use of a generalized Ito formula usually referred in the literature to as the Ito-Wentzell formula (see theorem 3.3.1, Kunita, 1990). This extended Ito formula incorporates in the same way as the classical Ito formula333relevant only to express the differential of a deterministic function of a stochastic process quadratic variation terms related to process but also co-variation terms between and the gradient of the random function . Its expression is in our case given by

 dϕ(t,Xt)= dtϕ+∇ϕ⋅dXt+∑id⟨\upartialϕ\upartialxi,(σ^Bt)i⟩dt+12∑i,jd⟨Xit,Xjt⟩\upartial2ϕ\upartialxi\upartialxjdt = 0, (21)

It can be immediately checked that for a deterministic function or for uncertainties directed along the function iso-values, the standard Ito formula is recovered. Otherwise, for a fixed grid point, function is hence solution of a stochastic differential equation of the form

 dtϕ(x,t)=v(x,t)dt+f(x,t)⋅d^Bt. (22)

The quadratic variation term involved in (21) is given through (14) as

 d⟨Xit,Xjt⟩=aij(x,t), (23)

 d⟨\upartialϕt\upartialxi,(σ^Bt)i⟩ =limn→∞1n∑jn∑ℓ=0(\upartialfj\upartialxi(x,∙,t)⋆φν(xℓ))σij(x,∙,t)⋆φν(xℓ), =∑j\upartialfj\upartialxi⋆φν(x,t)σijφν(x,t). (24)

From these expressions and identifying first the Brownian term and next the deterministic terms of equations (21) and (22), we infer

 f(x,y,t)T =−∇ϕ(x,t)Tσ(x,y,t), v(x,t) =−∇ϕ⋅w+∑i,j12aij\upartial2ϕ\upartialxi\upartialxj+∇ϕ⋅\upartialσ∙jν\upartialxiσijν,

and finally get

 dtϕ =Lϕdt−∇ϕ⋅σd^Bt, (25) Lϕ =−∇ϕ⋅w+∑i,j12aij\upartial2ϕ\upartialxi\upartialxj+∇ϕ⋅\upartialσj∙φν\upartialxiσijφν.

This differential at a fixed point, , defines the equivalent in the deterministic case of the material derivative of a function transported by the flow. For interested readers, several conditions under which linear equations with a multiplicative noise exhibits a unique weak, mild or strong solution are studied in (Prato and Zabczyk, 1992). In our case, only existence of a weak solution is needed as we will proceed to a formal integration by part. The differential of the integral over a material volume of the product is given by

 d∫V(t)qϕ(Xt,t)dx =d∫Ωqϕdx, =∫Ω(dtqϕ+qdtϕ+dt⟨q,ϕ⟩)dx,

where the first line comes from and the second one from the integration by part formula of two Ito processes. Hence, from (25) this differential is

 ∫Ω(dtqϕ+qLϕ+∇ϕ⋅a∇q)dtdx−∫Ωq∇ϕ⋅σd^Bt.

Introducing the (formal) adjoint of the operator in the space with Dirichlet boundary conditions, this expression can be written as

 ∫Ω((dtq+L∗q−∇⋅(a∇q))dt+∇⋅(qσd^Bt))ϕdx.

With the complete expression of and if , where stands for the characteristic function, we get the extended form of this differential:

 (26)

which can be further simplified as

 (27)

At last this expression may be more compactly written as

 ∫V(t)[dtq+∇⋅(qw)dt−12∑ij\upartial2\upartialxi\upartialxj(aijq)|∇⋅σ=0dt+12∥∇⋅σ∥2qdt+∇⋅(qσd^Bt)]dx,

where the third term must be computed by considering a divergence free diffusion tensor. Let us note it is now straightforward to get the expression of the transport theorem for a divergence free turbulent component:

 (28)

## 4 Mass conservation

This expression of volumetric rate of change relation allows us stating the mass conservation principle under fluid particles location uncertainty. Applying the previous transport theorem to the fluid density , we get a general expression of the mass variation:

 ∫V(t)[dtρ+(∇⋅(ρw)−12∑i,j\upartial2\upartialxi\upartialxj(aijρ)|∇⋅σ=0+12∥∇⋅σ∥2ρ)dt+∇⋅(ρσd^Bt)]dx. (29)

A mass conservation constraint on the transported volume implies canceling this expression. Besides, as the volume is arbitrary this provides the following constraint

 dtρ+∇⋅(ρw)dt=12(∑i,j\upartial2\upartialxi\upartialxj(aijρ)|∇⋅σ=0−12∥∇⋅σ∥2ρ)dt−∇⋅(ρσd^Bt). (30)

### 4.1 Incompressible fluids

For an incompressible fluid with constant density, canceling separately the slow deterministic terms and the highly oscillating stochastic terms of the mass conservation constraint (30), we obtain {multiequations}

 ∇⋅(σd^Bt)=0,∇⋅w=12∑i,j\upartial2aij\upartialxi\upartialxj. (31)

In the very same way as in any large scale flow decompositions, we aim here at recovering a physically plausible volume preserving drift component. Imposing to the drift component to be divergence free, this system boils down hence to an intuitive system of incompressibility constraints {multiequations}

 ∇⋅(σd^Bt)=0,∇⋅w=0, (32)

coupled with a less intuitive additional constraint on the quadratic variation tensor:

 ∇⋅(∇⋅a)=0. (33)

For the Kraichnan model (or for any divergence-free homogeneous random fields) this last constraint is naturally satisfied, as its quadratic variation is constant and the system comes down to the classical incompressibility constraint. In the same way, any homogeneous divergent-free random fields associated to a volume preserving drift leads to mass conservation.

### 4.2 Isochoric flows and isoneutral uncertainty

For divergence-free volume preserving flows (refered as isochoric flow) with varying density we get a mass conservation constraint of the form

 dtρ+∇ρ⋅wdt−12∑i,j\upartial2\upartialxixj(ρaij)dt=∇ρ⋅σd^Bt. (34)

An interesting property emerges if the uncertainty has a much larger scale along the density tangent plane than in the density gradient direction.This situation occurs in particular in oceanography or in meteorology where the fluid is stratified and the horizontal scale much larger than the vertical scale. In a such case, it is essential that the large-scale numerical simulations preserve the natural fluid stratification and consequently to define subgrid models with controlled diffusion along the iso-density surfaces. This behavior can be easily set up by defining the diffusion tensor as a diagonal projection operator so that the uncertainty lies on the isodensity surfaces:

 σij=δij−\upartialxiρ(x)\upartialxjρ(y)∥∇ρ∥2δ(x−y). (35)

Together with the small slope assumption Gent and McWilliams (1990) () the diffusion tensor and the quadratic variation can then be written as a matrix function

 a(x)=⎛⎜ ⎜⎝10αx(x)01αy(x)αx(x)αy(x)|α(x)|2⎞⎟ ⎟⎠, (36)

where we introduced the neutral slope vector . Owing to the divergence free condition of the diffusion tensor, this slope vector is necessarily constant along the depth axis and divergence free . This yields for the mass preserving equation a diffusion along the density tangent plane. The density evolves as the deterministic model:

 \upartialρ\upartialt+∇ρ⋅w=12∑ij\upartialxi(aij\upartialxjρ), (37)

where we considered the divergence free constraint . This type of anisotropic diffusion corresponds exactly to the so-called isoneutral diffusion currently used to model unresolved mesoscale eddies in coarse scale resolution of ocean dynamics simulations (Gent and McWilliams, 1990). As mentioned previously, the divergence free uncertainty induces a slope vector that is also divergence free and independent of the depth. This enforces further the density surfaces to be quasi-harmonic (for smooth depth density variation) in the horizontal planes and provides a strong stratification of the density organized as a pile of pancakes with no maxima excepted on the domain horizontal boundaries. To our knowledge those additional constraints are not taken into account in the isoneutral diffusion and according to our interpretation this comes to consider uncertainties that are not volume preserving.

In the case of the Kraichnan model the density variation is simplified as

 dtρ+∇ρ⋅wdt−γ12∑i,j∇2ρdt=∇ρ⋅σd^Bt, (38)

and when the drift term does not depend on random argument (i.e. if it corresponds to the flow expectation as in the case of a mean field dynamics) the density expectation evolution is a classical intuitive advection diffusion equation

 \upartialt¯ρ+∇¯ρ⋅w=12γ∇2¯ρ. (39)

The mass conservation constraint and the stochastic version of the Reynolds theorem allows us now expressing the linear momentum conservation equations.

## 5 Conservation of momentum

Newton’s second law states that the net force acting on the fluid is equal to the rate of change of the linear momentum. In our case, the flow evolution is described through an Ito diffusion (1), where the drift component represents the unknown flow we wish to estimate whereas the noise term denotes the fluctuating part either caused by physical processes that are neglected or generated by numerical errors and coarsening processes. Whatever the type of uncertainties considered on the fluid particles location, they introduce de facto an uncertainty on the forces. In order to take into account all the uncertainty sources within the momentum equations we consider a stochastic conservation principle of the form

 d∫V(t)ρ(w(x,t)dt+σ(x,t)d^Bt)dx=∫V(t)F(x,t)dx. (40)

In this equation the acceleration is highly irregular and has to be interpreted in the sense of distribution. For every :

 ∫h(t)∫V(t)F(x,t)dxdt=−∫h′(t)∫V(t)σ(x,t)d^Btdxdt+∫h(t)d∫V(t)ρw(x,t)dxdt. (41)

Since both sides of this equation must have the same structure, the forces can be written as

 ∫h(t)∫V(t)F(x,t)dxdt=−∫h′(t)∫V(t)σ(x,t)d^Btdxdt+∫h(t)∫V(t)(f(x,t)dxdt+θ(x,t)d^Bt)dx. (42)

The first terms of the right-hand sides of equations (41) and (42) are identical and cancel out. The second right-hand term of equation (41) corresponds to the derivative of the momentum associated to the resolved velocity component, whereas the second right-hand term of (42) provides us the structure of the forces under localization uncertainty. We now further develop those different terms.

According to the previous section results, the transport equation applied to the linear momentum gives for each component of the velocity:

 d∫V(t)ρwidx=∫V(t)[(dt(ρwi)+∇⋅(ρwiw)+∥∇⋅σ∥2ρwi−∑j,k12\upartial2\upartialxj\upartialxk(ajkρwi)|∇⋅σ=0)dt+∇⋅(ρwiσd^Bt)]dx. (43)

As for the forces, we will consider that only body forces and surface forces are involved (i.e., there is no external forces). A direct extension of the deterministic case provides us the surface forces as

 Σ=∫V−∇(pdt+d^pt)+μ(∇2U+13∇(∇⋅U)),

where is the dynamic viscosity, denotes the deterministic contribution of the pressure and is a zero mean stochastic process describing the pressure fluctuations due to the random velocity component. The gravity force is deterministic and standard.

Expressing the balance between those forces and the acceleration (43), incorporating then the mass preservation principle (30), and finally equating, in the one hand, the slow temporal bounded variation terms and, in the other hand, the Brownian terms, we get the expression of the flow dynamics:

 (\upartialw\upartialt+w∇Tw)ρ−12∑i,jaijρ\upartial2w\upartialxi\upartialxj− ∑i,j\upartial(aijρ)\upartialxj|∇⋅σ=0\upartialw\upartialxi =ρg−∇p+μ(∇2w+13∇(∇⋅w)), (44a) ∇d^pt=−ρw∇Tσd^Bt+μ(∇2 (44b) dtρ+(∇⋅(ρw)−12∑i,j\upartial2\upartialxi\upartialxj(aijρ) |∇⋅σ=0+12∥∇⋅σ∥2ρ)dt=∇⋅(ρσd^Bt). (44c)

This system provides a general form of the Navier Stokes equations under location uncertainty. Similarly to the Reynolds decomposition, the dynamics associated to the drift component includes an additional stress term that depends on this resolved component and on the uncertainty diffusion tensor. In order to specify further the different terms involved in this general model, we will examine simpler particular instances of this system. In the following, we will confine ourselves to the case of an incompressible fluid of constant density. Considering as a first example the divergence free isotropic Kraichnan model defined in (17) and its associated transport equation, we obtain straightforwardly the following Navier-Stokes equation:

 (\upartialw\upartialt+w∇Tw−γ12∇2w)ρ=ρg−∇p+μ∇2w, (45a) ∇d^pt=−ρ(w∇T)dξt+μ∇2dξt, (45b) ∇⋅w=0. (45c)

The first equation of this model corresponds to a traditional turbulent diffusion modeling relying on the Boussinesq assumption with a constant eddy viscosity coefficient. Let us note that for this system the second equation related to the random pressure term is not needed to compute the resolved drift . The uncertainty is in that case specified a priori as an homogeneous random field. Considering a more general divergence free random component, we obtain {multiequations}

 \singleequation(\upartialw\upartialt+w∇Tw−12∑i,j\upartialxi\upartialxj(aijw))ρ=ρg−∇p+μ∇2w, (46)
 \singleequation∇d^pt=−ρ(w∇T)σd^Bt+μ∇2σd^Bt, (47)
 \tripleequation∇⋅w=0,∇⋅σ=0,∇⋅(∇⋅a)=0. (48)

This model of Navier-Stokes equations includes now a diffusion term that cannot be directly related to the traditional Boussinesq eddy viscosity formulation anymore. Furthermore, the incompressibility condition on the variance tensor provides a non-local constraint on the subgrid model, which is lacking in usual eddy viscosity models (Kraichnan, 1987). One can point out that for divergence-free incompressible uncertainty models this term is globally dissipative. As a matter of fact the energy associated to the subgrid term reads

 ∫ΩwT∑i,j\upartial2\upartialxi\upartialxj(aijw)= −∫Ω∑k∇wTka∇wk−∫Ω(∇⋅a)∑k∇wkwk = −∫Ω∥∥∇w∥∥2a−∫Ω(∇⋅a)∇(12∥w∥2). (49)

The first term , is always non-negative as, , is semi-positive definite. The second term associated to the incompressibility constraint, , cancels:

 ∫Ω(∇⋅a)∇(12∥w∥2)dx =12∫Ω(∇⋅(∇⋅a∥w∥2))dx (50) =12∫\upartialΩ∥w∥2(∇⋅a)⋅nds=0.

The energy of the subgrid term reduces to