# Stochastic representation of the Reynolds transport theorem: revisiting large-scale modeling

## Abstract

We explore the potential of a formulation of the Navier-Stokes equations incorporating a random description of the small-scale velocity component. This model, established from a version of the Reynolds transport theorem adapted to a stochastic representation of the flow, gives rise to a large-scale description of the flow dynamics in which emerges an anisotropic subgrid tensor, reminiscent to the Reynolds stress tensor, together with a drift correction due to an inhomogeneous turbulence. The corresponding subgrid model, which depends on the small scales velocity variance, generalizes the Boussinesq eddy viscosity assumption. However, it is not anymore obtained from an analogy with molecular dissipation but ensues rigorously from the random modeling of the flow. This principle allows us to propose several subgrid models defined directly on the resolved flow component. We assess and compare numerically those models on a standard Green-Taylor vortex flow at Reynolds 1600. The numerical simulations, carried out with an accurate divergence-free scheme, outperform classical large-eddies formulations and provides a simple demonstration of the pertinence of the proposed large-scale modeling.

## 1Introduction

The large-scale modeling of fluid flow dynamics remains nowadays a major research issue in fluid mechanics or in geophysics despite an enormous research effort since the first investigations on the subject 150 years ago [5]. The research themes behind this topic cover fundamental issues such as turbulence modeling and the analysis of fully developed turbulent flows, but also more applicative research problems related to the definition of practical numerical methods for the simulation of complex flows. In this latter case the difficulty consists in setting up a reliable modeling of the large-scale dynamics in which the contribution of unresolved small-scale processes are explicitly taken into account. For the Navier-Stokes equations, the problem is all the more difficult that the spatial and temporal scales are tightly interacting together.

The neglected processes include, among others things, the action of the unresolved motion scales, complex partially-known forcing, an incomplete knowledge of the boundary conditions and eventual numerical artifacts. Such unresolved processes must be properly taken into account to describe accurately the energy transfers and to construct stable numerical simulations. In real world situations, the complexity of the involved phenomenon prevents the use of an accurate – but inescapably expensive – deterministic modeling. We advocate instead the use of a stochastic modeling.

Within this prospect, we aim at describing the missing contributions through random fields encoding a flow component only in a probability distribution sense. Those variables correspond to the discrepancies or errors between the dynamical model and the actual dynamics. Their modeling is of the utmost importance in geophysics, either for data assimilation or forecasting issues. In both cases, an accurate modeling of the flow errors dynamics enables to maintain an ensemble of flow configurations with a sufficient but also meaningful diversity.

Small-scale processes are responsible both for an energy dissipation but also for local backscattering of energy [53]. The introduction of random variables in the flow dynamics has been considered by several authors, as it constitutes an appealing mechanism for the phenomenological modeling of intermittent processes associated to the inverse energy cascade [36]. Recently those models have regained a great interest for the modeling of geophysical flows [9] in climate sciences (see also the thematic issue [50] or the review [16]).

Numerous turbulence models proposed in the context of Large Eddies Simulations (LES) and Reynolds Average Simulations (RANS) introduce *de facto* an eddy viscosity assumption to model the energy dissipation due to unresolved scales [14]. This concept dates back to the work of Boussinesq [5] and Prandtl [56]. It relies on the hypothesis that the energy transfer from the resolved scales to the subgrid scales can be described in a similar way to the molecular viscosity mechanism. It is therefore not at all related to any uncertainty or error quantities. In models dealing explicitly with a statistical modeling of the turbulent fluctuations there is thus some incoherency in representing directly the dissipative mechanism attached to random terms through an eddy viscosity assumption. In this work we will not make use of such an hypothesis. Instead, we will rely on a general diffusion expression that emerges naturally from our formalism.

This subgrid model is properly derived from a general Lagrangian stochastic model of the fluid motion in which the fluid parcels displacement is decomposed in two components: a smooth differentiable (possibly random) function and a random field, uncorrelated in time but correlated in space. Such a decomposition consists in separating or “filtering” a rough velocity in a smooth slow time-scale component and a fast oscillating velocity field representing the unresolved flow. Though there is, in general, no sharp time-scale separation in turbulent flows, the resolved velocity can be interpreted as a temporally coarse-grained component whereas the time-uncorrelated component stands for the small time-scale unresolved velocity. As a temporal smoothing imposes implicitly a spacial smoothing, this separation can be thus interpreted in terms of a LES filtering technique. Yet, the corresponding Eulerian formulation does not ensue from a filtering procedure. It is thus not prone to errors associated to the violation of the commutation assumption between the filter and the spatial derivatives [21]. Besides, those equations introduce an effective advection related to the small-scale velocity inhomogeneity. This modified advection, empirically introduced in Langevin models of particle dispersion [39], corresponds exactly to a phenomenon, termed *turbophoresis*, related to the migration of inertial particles in regions of lower turbulent diffusivity [64].

The large-scale representation of the Navier-Stokes equations on which we rely in this study are built from a stochastic version of the Reynolds transport theorem [43]. This modified Reynolds transport theorem, which constitutes here the cornerstone of our large-scale dynamics representation, is presented in the following section. General invariance properties of the corresponding large-scale dynamics such as scale and Galilean invariances are detailed in a comprehensive apppendix. In Section 3 several novel subgrid tensors will be devised and compared on a standard Green-Taylor vortex flow [68]. We will show that all the proposed schemes outperform the usual dynamic Smagorinsky subgrid formulation [19].

## 2Stochastic modeling of fluid flow dynamics

Numerous methodological choices can be envisaged to devise stochastic representations of the Navier-Stokes equations. The simplest method considers additional random forcing to the dynamics. This is the choice that has been the most often performed since the work of Benssoussan [3]. Another choice, in the wake of Kraichnan’s work [29], consists in closing the large-scale flow representation in the Fourier space by relying on a Langevin equation [31]. Obviously the frontiers between those two methodologies are sometimes fuzzy, and numerous works rely on both strategies in order to devise the shape that should take the random variables evolution law [25]. Lagrangian stochastic models based on Langevin equations have been also intensively used for turbulent dispersion [62] or in probability density function (PDF) modeling of turbulent flows [24]. Those Lagrangian models, which require to model the drift and diffusion functions, lead to very attractive particle based representations of complex flows [54]. They are nevertheless not adapted to global large-scale Eulerian representations of the flow dynamics.

In this work, we will rely on a different framework in specifying the stochastic nature of the velocity from the very beginning as proposed in [43]. The basic idea is built on the assumption that the Lagrangian fluid particles displacement results from a smooth velocity component and a highly oscillating stochastic velocity component uncorrelated in time,

with the velocity components:

In this decomposition the first right-hand term is a smooth function of time associated to the large-scale velocity component. The second term stands for the small-scale velocity field. It is a white noise velocity component defined from the (formal) time-derivative of the random field: . This random field is a three-dimensional centered Wiener process; it is thus uncorrelated in time but can be anisotropic and inhomogeneous in space. Since we focus in this study only on incompressible flows, the small-scale component is defined as a divergence-free random field; it is hence associated to a divergence-free diffusion tensor:

Analogously to the standard deterministic case, the derivation procedure from the physical conservation laws of the Navier-Stokes equations is based primarily on the Reynolds transport theorem (RTT).

### 2.1Stochastic Reynolds transport theorem

The RTT provides the expression of the rate of change of a scalar function, , within a material volume, . For a stochastic flow (Equation 2) with an incompressible small-scale velocity component (), this expression derived in [43] (see also Appendix A where a comprehensive derivation is provided for readers convenience), is given by:

This modified RTT involves the time increment of the random scalar quantity (the differential of at a fixed point) instead of the time derivative. A diffusion operator emerges also naturally. For clarity’s sake, this term is designated as “subgrid stress tensor” following the protocols of large eddies simulation (LES). However, its construction is quite different. It is not based on Boussinesq’s eddy viscosity assumption nor on any structural turbulence models [61] but arises directly from stochastic calculus rules. It expresses the mixing process exerted on the scalar quantity by the fast oscillating velocity component. This diffusion term is directly related to the small-scale component through the *variance tensor*, , defined from the diagonal of the small-scale velocity covariance:

it can be checked that the variance tensor corresponds to an eddy viscosity term (with units in . This term plays thus a role similar to the eddy viscosity models introduced in classical large scale representations [2] or to the spectral vanishing viscosity [28]. It is also akin to numerical regularization models considered in implicit models [1]. Our approach is nevertheless more general as it does not rely on *a priori* fixed shape of the Reynold stress (e.g. Boussinesq assumption) nor does it presuppose a given numerical discrete scheme (e.g. implicit models).

A corrective advection term, , appears also in the stochastic RTT formulation. This correction expresses the influence of the small scales inhomogeneity on the large ones. A drift is engendered from the regions associated with maximal variance (maximal turbulent kinetic energy - TKE) toward the area of minimum variance (e.g. minimal TKE). This corresponds to the *turbophoresis* phenomenon associated with turbulence inhomogeneity, which drives inertial particles toward the regions of lower diffusivity [7]. For homogeneous noise, the variance tensor is constant and this corrective advection does not come into play. It can be noted that this advection correction is of the same form as the one proposed in [39].

Through this modified RTT, stochastic versions of the mass and momentum conservation equations can be (almost) directly derived. Incompressibility conditions can for instance be immediately deduced from the RTT applied to and the flow jacobian ():

Together with the incompressibility of the random term, the incompressibility condition reads thus:

In the case of an incompressible large-scale flow component, , this reduces to:

Note that for a divergence-free isotropic random field such as the Kraichnan model [30] the last condition is naturally satisfied, since this unresolved velocity component is associated with a constant variance tensor.

### 2.2Mass conservation

Applying the RTT to the fluid density, , and canceling this expression for arbitrary volumes, we get the following mass conservation constraint:

For an incompressible fluid with constant density, together with a volume-preserving (isochoric) condition on the large-scale velocity component, we retrieve the incompressibility conditions (Equation 7). It can be noted also that equation (Equation 8) still constitutes a transport equation since it preserves energy. As a matter of fact, it can be shown that the energy intake brought by the small-scale component is exactly compensated by the energy dissipation associated to the diffusion term [58].

### 2.3Linear momentum conservation

The application of the stochastic version of the RTT on the stochastic momentum and the introduction of the forces acting on the flow enables to derive from the second Newton law the following Navier-Stokes equations [43]:

In this expression is the dynamic viscosity, denotes the large-scale (slow) pressure contribution and is a zero-mean turbulent pressure related to the small-scale velocity component. Similarly to the classical Reynolds decomposition, the dynamics of the resolved component includes an additional stress term that depends here on the resolved velocity component. A correction of the advection velocity also occurs. Both terms depend on the variance tensor which gathers the action of the turbulent fluctuations on the large-scale velocity.

It can be observed that the large-scale energy evolution is dissipative. This generalizes thus the Boussinesq 1877 assumption, which conjectures a dissipative effect of the turbulent fluctuations on the mean flow. In the case of a divergence-free isotropic random field (with a constant diagonal variance tensor), this system boils down to an intuitive constant eddy viscosity diffusivity model:

where the Laplacian dissipation term is augmented by the random field variance. The use of constant eddy viscosity thus finds here its justification as a direct consequence of an isotropic turbulence assumption.

The subgrid stress tensor involved in our formalism constitutes an anisotropic diffusion whose preferential directions of diffusion are given by the small-scale velocity variance tensor. Setting the diffusion tensor, , or its associated variance tensor allows us to define directly the subgrid diffusion term and the effective advection. For instance imposing to the small-scale random velocity to live on the iso-density surfaces provides immediately a clear justification of the isopycnal diffusion employed in oceanic circulation models [43]. The specification of the turbulent fluctuations in terms of a stochastic process provides a means to interpret different subgrid models but also to devise new ones either through *a priori* specifications or data-driven strategies.

General invariance properties of the proposed large-scale representation are listed in the Appendix C. We briefly summarize them here. It is in particular shown that the distribution of the velocity anomaly is in the general case not Gaussian and does not consequently correspond to normal or quasi-normal approximations [46]. We show that this stochastic representation has remarkable invariance properties; it is Galilean invariant and preserves (in the absence of molecular viscosity) the Euler equations’ scale invariance properties. Otherwise, a useful scaling for the variance tensor is derived from the Kolmogorov-Richardson scaling and a dimensionless number relating the large-scale kinetic energy to a characteristic value of the velocity variance at the resolved scale.

Interested readers may also consult [58] for the derivation of stochastic representation of several geophysical flows.

## 3Numerical simulation and assessment

In this section we assessed, through numerical simulations, the performances of the proposed large-scale dynamics for different variance tensor models. Those simple models have been defined from local statistics of the resolved component and compared to the classical Smagorinsky subgrid model associated with a dynamical procedure [19]. Those numerical experiments have been performed on the Taylor-Green flow [68].

### 3.1Taylor-Green vortex flow simulation

Taylor-Green vortex flow is a critical test for numerical schemes, as both the convective term and viscous term play important roles. Due to the energy cascade generated by the convective term, the flow becomes rapidly turbulent with the creation of small-scale structures up to a dissipation peak. This stage is followed by a decay phase similar to a decaying homogeneous turbulence. As a consequence, a precise and high-order representation of the viscous and convective terms is needed to get an accurate numerical simulation. This flow is considered as a prototypical system to study the production of small-scale eddies due to vorticity increase and vortex stretching mechanism [6].

In Cartesian coordinates, this flow is defined by the following initial conditions:

and

The computation domain is defined as a cubic box with periodic boundary conditions on all the faces. The length of the domain is set to in each of the axis direction, which gives a characteristic length and a Reynolds number . In the literature, several high-order numerical methods have been tested for the direct numerical simulation (DNS) of the Taylor-Green vortex flow, see [11] and references therein. In this study we used a discrete scheme built on a divergence-free wavelet basis [26]. This scheme presents several computational advantages. First of all, it guarantees a divergence-free solution in the physical domain with a good numerical complexity. Besides, as the spatial filters considered corresponds to a multi-resolution projection, two successive filtering operations can be switched together. This property reveals useful within the Germano dynamic strategy enabling to estimate the subgrid tensor weight factor. This numerical scheme achieves similar performances to a pseudo-spectral method.

### 3.2Analysis criterion

The different numerical simulations performed are mainly assessed and compared according to the evolution along time of the following criterion:

The mean kinetic energy (KE)

The mean kinetic energy dissipation rate

where is the rate of strain tensor: . The mean kinetic energy is linked to the dissipation rate by:

To clearly separate those two forms of the energy dissipation we will denote the expression computed from the kinetic energy differentiation and the dissipation computed from the rate of strain norm.

Besides, in the sequel, all those averaged quantities are computed in their dimensionless form:

where denotes the convective time. The temporal evolution of those mean energy quantities enable to monitor the quality of the solution over time and to assess the accuracy of the discrete scheme used for the velocity gradients. In addition to those criterion, we will plot the energy spectrum of the resolved velocity at several distinct times.

### 3.3Direct Numerical Simulation of the Taylor-Green flow

We evaluated first the ability of the divergence-free wavelet based method to reproduce faithfully the main features of the Taylor-Green flow. For this purpose, two different direct numerical simulations have been conducted.

The first experiment concerns a simulation at a low Reynolds number: . For this case the divergence-free wavelet has been run on a regular grid of points. For comparison purpose we performed a classical pseudo-spectral simulation with Fourier modes. Let us note that at this low Reynolds number (), only Fourier modes are required to represent accurately all the hydrodynamics scales in the limit of , where is the Kolmogorov scale [11]. The time evolution of the mean kinetic energy and the mean dissipation rate obtained for both methods are plotted on the left and right panels of Figure 1 respectively. As can be observed on those two graphics the solutions superimpose almost perfectly. For both methods the energy dissipation computed from the rate of strain tensor norm, , and from the kinetic energy differentiation, , coincide perfectly. We plotted therefore only the rate of strain norm.

In a second experiment, both simulations have been then carried out for a moderate Reynolds number fixed at . Since the study of [6], numerous numerical experiences have been conducted to understand accurately the dynamic of the Taylor-Green vortex flow at this Reynolds number. As mentioned in [11], all the scales of the flow are captured with degrees of freedom but degrees of freedom are sufficient to represent its main characteristics. We therefore run the divergence-free wavelet based simulation on a mesh grid. We also performed a Fourier pseudo-spectral simulation with a dealiasing procedure at the same resolution. On Figure 2 we displayed for both simulations the kinetic energy and dissipation rate time evolutions. On the same figure, for comparison purpose, we also plotted the curves corresponding to a Fourier modes pseudo-spectral solution (with dealiasing) available from [11]. As in the previous case the dissipation matches . We therefore only show the rate of strain tensor norm in the right panel of Figure 2.

It can be observed that the divergence-free wavelet simulation is in good agreement with the spectral reference solution, especially before the dissipation peak when the convective phenomena is predominant. Slight discrepancies appear in the decaying phase starting after the dissipation peak ().

Those experiments show that the divergence-free wavelet based method provides results of comparable quality to those obtained from pseudo-spectral simulations at identical resolutions. The wavelet simulation performed on the mesh grid will serve as a reference “DNS” for comparison purpose. We relied on this method to carry out all the numerical simulation performed in this study. The performances of various subgrid modeling are discussed in the next section.

### 3.4Variance tensor and subgrid modelling

One of the main advantages of the stochastic formalism we propose lies in the great flexibility of the anisotropic diffusion specification. The variance tensor, , can be fixed from *a priori* shape imposed either directly on the small scale variance or on the diffusion tensor. In some cases, such knowledge could probably be inferred from the physical approximations considered to constitute the model. Aspect ratio simplifications and/or boundary conditions that are not perfectly known could be used as well to constrain the small-scale velocities to live on specific iso-surface. Another route would consist in specifying the turbulent velocity components from small-scale measurements statistics. Despite all those directions are worth exploring as they open new strategies to design sub-grid tensors, in this study we choose to focus on several simple models of tensor to explore the potential of the proposed modeling.

#### Empirical specification through local mean

The first model consists in assuming ergodic assumption to compute the variance tensor from local statistics of the resolved velocity component. The variance tensor is here defined from empirical velocity covariance computed on a local neighborhood:

where denotes the empirical mean on a spatial or temporal window , is the spatial scale considered for the simulation, which corresponds to a mesh grid composed of grid points, and is the Kolmogorov scale. The empirical averaging, , is computed either spatially over a small () window centered around point or temporally, at point , over the time interval . In the following, they are referred to as the spatial and temporal local covariances respectively. Both models are weighted by the variance tensor scale ratio derived in section ?.

#### Optimal specification through scale similarity

The second model is defined from two successive filtering of the resolved velocity component and a scale similarity assumption. The filtering is defined through the associated wavelet multi-scale projector. More precisely, since the velocity can be decomposed as:

with the projector onto the scaling functions basis [15], the resolved filtered velocity is defined by:

The second, so called “test”, filtering is here defined by the projection of onto resolution , which is the immediate coarser resolution following the simulation scale:

The stochastic Navier-Stokes equations ( ?) rewritten respectively for and reads:

and

Taking the difference of the momentum equations and for the two subsequent levels and provides the residual dynamics:

where , and . Note that a similarity assumption has been imposed for the variance tensor. It can be noticed that, if the projector commutes with differentiation, due to the filtering projection property, the Stokes equation in the left-hand side cancels after filtering. Then, instead of solving the right-hand expression at level and then projecting back at a finer level the estimated variance tensor, we rather propose to solve it at scale, , in a least squares sense. Introducing the variance tensor incompressibility constraint we seek the minimizer of the following nonlinear functional:

with the positivity constraint:

In practice the minimization of the functional has been carried out using a quasi-newton method combined with the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method [47] to approximate the Hessian matrix. The variance tensor has been assumed to be diagonal at all points. In order to impose the positivity constraint, instead of computing we preferred to compute its square root . For details on the computation of the gradient of , we refer to Appendix D.

#### Smagorinsky subgrid model

The third model evaluated corresponds to the classical Smagorinsky eddy viscosity formulation,

coupled with the Germano estimation procedure [19] to fix dynamically the eddy viscosity constant from a least squares estimation [38] and a filtering of the velocity field at two consecutive scales. The “test” filtering is as previously defined by the projection of onto the immediate coarser resolution . Denoting the filtering operation , the constant is given by:

with

and denotes a spatial averaging.

### 3.5Large-scale simulation numerical results

The performances of the different sub-grid models have been compared for a mesh grid, which corresponds to a resolution scale . Figure 3 shows for the different models the time evolution of the mean kinetic energy (left panel) and the mean energy dissipation rate (right panel). Those curves can be compared to the reference solution computed on a grid (). An important remark can be raised here. The filtered velocity fields do not correspond to the solution of a Navier-Stokes equation started from the large-scale (filtered) initial conditions. It is thus not relevant to correlate the large-scale solutions with the filtered velocity. It is more meaningful to compare them directly with the DNS.

From the rate of strain tensor norm displayed on the right side of Figure 3, two groups of subgrid tensor can be recognized. The first one gathers the dynamic Smagorinsky tensor (`Dyn-S`

) and the variance tensor defined from the spatial local velocity covariance (`Spatial`

). The second group is composed of the optimal variance tensor with a scale similarity assumption (`Opt-SS`

) and the variance tensor built from the velocity temporal covariance (`temporal`

). For both groups, we displayed on the right-hand side of Figure 4 and Figure 5 the energy dissipation rate computed from the kinetic energy rate of change. On the left-hand side we plotted the kinetic energy associated to each model. Those curves can be compared to the DNS reference curve. We observe from those figures that all the large-scale models achieved to reproduce the right amount of kinetic energy and energy dissipation rate. The dynamic Smagorinsky subgrid tensor is the model which provides the lower kinetic energy and the lower rate of strain norm. We note the first group of subgrid tensors exhibits in a general way a too fast dissipation in the first phase of the flow (), characterized by the domination of the advection mechanism. An undue smoothing of the resolved velocity gradient results from this too high dissipation rate. This is confirmed by looking at the rate of strain norm plotted in Figure 3. Within this first group, the model build from local spatial velocity covariance performs better than the original Smagorinsky tensor. The second group of variance tensor outperforms clearly the results of the first group. Both the temporal covariance model and the optimal variance tensor show strikingly very close results. The optimal variance tensor has a slightly higher dissipation rate at the dissipation peak whereas the temporal covariance model fits almost perfectly the DNS results around the dissipation peak (). Both models slightly differ from the DNS in the decaying phase. In a general way, we note that all the models have a too slow dissipation rate at the end of the decaying phase (). They all smooth too much the velocity gradients in that regime. It must be outlined that the good behavior of the temporal covariance confirms the relevance of the variance tensor scaling since no dynamics strategy has been performed in that case.

We display next on Figure 6 and on Figure 7 the energy spectrum associated to the different models at four different instants. All the models provide satisfying solutions with similar spectrum. Compared to the other models, the dynamic Smagorinsky subgrid stress produces a noticeable energy bump at the cutoff scale. The temporal and optimal variance tensor have spectrum which are in general closer to the DNS spectrum. At the end of the turbulence decay phase () all the models provide close results.

As we rely on a wavelet scheme for the numerical simulations, it is insightful to inspect the discrete power spectra computed from the wavelet coefficients. They are plotted on figures Figure 8 and Figure 9. Wavelet power spectrum corresponds to an averaged version of the Fourier spectrum [8] and exhibits the same slopes as the Fourier spectrum [52]. A discrete version of the wavelet spectrum as plotted here provides one energy measure per scale level. It can be observed on the four spectra that the dynamic Smagorinsky tensor exhibits an undue energy intake at the cutoff scale. This amplification of energy is likely due to noisy velocity fields at the cutoff. The three different model ensuing from our statistical representation of the small-scale velocity component does not show such a deficiency.

For information purpose we draw on Figure 10 the value along time of the constant weighting the dynamic Smagorinsky model. The obtained maximum value is about and the mean value is , this is in good agreement with the predicted values, see [19].

Regarding the different criteria explored, we see that the proposed large-scale subgrid tensors based on local velocity covariances perform better than the dynamic Smagorinsky subgrid model. Among all the models proposed, the optimal variance tensor estimation based on a two-scale similarity assumption and the variance tensor constructed from a local temporal variance of the resolved velocity field outperform significantly the others. They lead to very good results that are strikingly close to each other. To our knowledge they both outperform the state-of-the-art solutions for this flow [11]. Though very simple, the temporal covariance model weighted by a unique constant fixed through the quadratic scaling rule presented section ? gives very good results. This demonstrates the pertinence of this scaling. The choices operated here are among the simplest that can be devised. More involved schemes could be easily imagined. For instance, variance tensors based on vorticity statistics might be very interesting to explore. Another route would be to elaborate this tensor from statistics extracted from measurements or DNS data. Large-eddies simulation models could likely be proposed in this spirit for non-homogeneous turbulence or boundary layer flows. Note also that both correlation schemes rely on a constant fixed through a rough dimensional scaling. A dynamical procedure in the same vein as the one used for the dynamics Smagorinsky subgrid model could be beneficial to get a finer estimate of the constant. As our purpose was here to bring a simple demonstration of the wide potentiality offered by the proposed stochastic modeling, we leave such potential improvements to future works.

## 4Conclusion

In this paper we have described a decomposition of the Navier-Stokes equations in terms of a temporally smooth velocity component and a fast oscillating random field associated to the unresolved flow component. This decomposition leads to a new large-scale representation paradigm, which can be interpreted as a large eddies simulation formalized through a time-scale separation. An advection correction and a subgrid diffusion term both emerge from this formalism. They encode respectively *turbophoresis* phenomenon and anisotropic mixing effect due to turbulence. The corresponding subgrid tensor enables generalizing the Boussinesq eddy viscosity concept on solid a theoretical ground. Such large-scale representation have been assessed on a Taylor-Green vortex flow. We compared different models of the variance tensor built from local averaging and a scale similarity least squares estimation procedure. The different numerical simulations outperform a standard dynamic Smagorinsky model based on Boussinesq’s eddy viscosity assumption. We believe those first results constitute are very encouraging as more advanced models could yet yield substantial improvements. As a matter of fact this formalism, built from physical conservation laws, paves the way to new possibilities to design efficient subgrid schemes. One could for instance explore models where the diffusion tensor is learned from DNS data or from small-scale observations. Another route of investigation consists in defining adapted basis for the small-scale random field from the fluctuations observed on two consecutive scales of the resolved tensor.

## AStochastic Reynolds transport theorem

We derive in this appendix the expression of the Reynolds transport theorem for the drift plus noise decomposition: In a Lagrangian stochastic picture, the infinitesimal displacement associated to a trajectory of a particle is noted:

The random field involved in this equation is defined for all points of the fluid domain, through the kernel associated to the diffusion operator

This operator is assumed to be of finite norm for any orthonormal basis of the associated Hilbert space ) and to have a null boundary condition on the domain frontier:

Function denotes a -dimentional Brownian function (see [57] for more information on infinite dimensional Wiener process) and the resulting -dimensional random field, , is a centered vectorial Gaussian function correlated in space and uncorrelated in time with covariance tensor:

We observe that those random fields have a (mean) bounded norm: , where the trace of the covariance tensor is given for any basis of as . In the following we will note the diagonal of the covariance tensor as: . Tensor, , that will be referred to as the variance tensor

is by definition a symmetric positive definite matrix at all spacial points, . This quantity corresponds to the time derivative of the so-called quadratic variation process:

The notation stands for the quadratic cross-variation process of and . This central object of Stochastic Calculus, can be interpreted as the covariance along time of the increments of and . The quadratic variation process is briefly presented in Appendix B.

The drift term, , of Lagrangian expression (Equation 13), represents the “smooth” – differentiable – part of the flow whereas the random part,

figures the small-scale velocity component associated to a much thinner time-scale. This component is modeled as a delta-correlated random field in time as it represents a highly non regular process at the resolved time scale. In this work, we assume that this small-scale random component is incompressible and therefore associated to a divergence free diffusion tensor: . This assumption, which obviously does not prevent the drift component (and therefore the whole field) to be compressible, leads to much simpler modeling and remains realistic for the models considered in this study.

Let us consider now a spatially regular enough scalar function of compact support, transported by the stochastic flow (Equation 13) and that vanishes outside volume and on its boundary . As this function is assumed to be transported by the stochastic flow, it constitutes a stochastic process defined from its initial time value :

We will assume that both functions have bounded spatial gradients. Besides, the initial function vanishes outside the initial volume and on its boundary. Let us point out that in this construction, 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

which here cancels as is transported by the flow. A separation between the slow deterministic terms and the fast Brownian term would yield a specific noise term orthogonal to . This would boil down to the deterministic case, 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 [32]. This extended Ito formula incorporates in the same way as the classical Ito formula a 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

It can be immediately checked that for a deterministic function, the standard Ito formula is recovered since the co-variations terms between and cancel.

It follows from Equation 14 that for a fixed grid point, function is solution of a stochastic differential equation of the form

where the Brownian term must compensate the Brownian term of (Equation 14). The quadratic variation term involved in (Equation 14) is given through (Appendix B) as

and the covariation term reads

In these expressions (resp. ) designates the kernel associated to operator (resp. ). Now, identifying first the Brownian term and next the deterministic terms of equations (Equation 14) and (Equation 15), we infer

and finally get

This differential at a fixed point, , defines the equivalent in the deterministic case of the material derivative of a function transported by the flow. The differential of the integral over a material volume of the product is given by

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

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

With the complete expression of (remarking that the second right-hand term of ? is self-adjoint) and if , where stands for the characteristic function, we get the sought form of this differential:

## BQuadratic variation and covariation

We recall here the notion of quadratic variation and co-variation, which play a central role in stochastic calculus. We will here restrict ourselves to standard Ito processes. Quadratic variation and co-variation correspond respectively to the variance and covariance of the process 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 :

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.

## CGeneral properties of the large-scale stochastic Navier-Stokes model

First of all, it is important to emphasize that the distribution of the velocity anomaly, , with the Eulerian velocity field

is, in the general case, *not* Gaussian. As a matter of fact, due to the multiplicative noise (the diffusion tensor depends *a priori* on the flow) and as the resolved dynamics is nonlinear, the ensemble mean of the flow velocity is in general not given by . Hence, this construction does not correspond to a turbulence model with a Gaussian closure hypothesis of the fourth-order correlation functions as considered in the Millionschikov hypothesis [46]. Nor is it based on quasi-normal approximations such as the EDQNM approaches [48], which opposite to the previous hypothesis ensures the positivity of the energy spectrum. In the proposed decomposition, the velocity is clearly a Markovian stochastic process, but its distribution is not Gaussian. Note also that in our approach no isotropic assumption on the increments nor on the random component are considered to define the subgrid tensor.

As already mentioned, the proposed stochastic representation can be interpreted as a temporal decomposition of the original Navier-Stokes equations. In the following sections we list several properties of this model. We begin first by a useful scaling relation of the variance tensor.

### c.1Variance tensor scaling

In the conditions of Kolmogorov-Richardson scaling, at scale , the velocity increments and the eddies turn-over time scale as and respectively, with denoting a constant energy dissipation rate across the inertial scales range. The turn-over time ratio for two different scales in this range, , exhibits a direct relation between a change of time scale and a change of spacial resolution. A coarsening in time yields thus a space dilation.

Besides, at scale , we have:

At the smallest scale, this quantity corresponds to the variance tensor . We have thus:

Within a given cell, , at scale , the energy of the small-scale random field is given as:

which from (Equation 17) gives us:

This relation provides us an expression of the scaling between the subgrid variance tensor at scale and the resolved velocity anomalies. It will serve us to impose a proper tuning of the subgrid term when the variance tensor is defined from the resolved velocity.

### c.2Adimensionalization of large-scale stochastic Navier-Stokes equations

Considering the scaled coordinates , , with velocity , pressure and a variance tensor , we have:

Assuming, as in the previous section, that the characteristic value of the resolved variance tensor, , depends linearly on the variance tensor at the smallest scale, , the dimensionless form of the large-scale stochastic Navier-Stokes equations reads:

where we introduced a dimensionless number

which relates the typical kinetic energy at scale to a characteristic value of the velocity variance, , at the same level. In this expression we considered from (Equation 18) a characteristic value of the resolved variance tensor

where we introduced a scale factor ratio, (with K(1)=1). As previously inferred, under the Kolmogorov scaling hypothesis this ratio scales as . Let us finally observe that for a typical variance tensor, , tending to zero, (Equation 19) tends to the original Navier-Stokes equations. It is the situation occurring when the variance tensor is negligible in front of the drift energy.

In the next section we will see that the large-scale model proposed conserves the invariance properties of the original Navier-Stokes equations.

### c.3Invariance of large-scale stochastic Navier-Stokes equations

Let us consider in the following the classical symmetry groups of the incompressible Navier-Stokes equations in a periodic domain and in the absence of external forcing.

#### Translation-invariance

Space translation invariance is achieved if the whole velocity field is still a solution of the Navier-Stokes equations. In our setting we have

Separating the Brownian component from the smooth term yields: , and hence , , , , , and translational invariance follows immediately.

#### Time-shift invariance

Time-shift invariance is obtained recalling that Brownian motion has itself a time-shift invariance property. This property allows us to write (in law) for :

which leads straightforwardly to time-shift invariance.

### c.4Rotational and reflectional invariance

These two symmetry groups correspond to a constant rotation or to a reflection of the coordinates system. The transformed coordinates are and . Note that , which yields and hence . The invariance of the transformed Navier-Stokes equations follows immediately:

#### Galilean transformations invariance

A model is said to be Galilean invariant if it is invariant with respect to an inertial transformation of the representation frames: . The velocity in this translated non-accelerating system of reference reads

with

From those equations we obtain:

As the pressure and viscous force are Galilean invariant and since , this shows the system is Galillean invariant.

#### Time reversal invariance

As with the original Navier-Stokes equation, the stochastic model is not invariant to a time reversal transformation: . This property is respected only for the Euler equation. Despite the time reversibility property of Brownian motion this loss of symmetry is due to the even nature of the quadratic tensor.

#### Scale invariance

Considering the scaled transformation:

we get: