Stochastic representation of the Reynolds transport theorem: revisiting largescale modeling
Abstract
We explore the potential of a formulation of the NavierStokes equations incorporating a random description of the smallscale 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 largescale 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 GreenTaylor vortex flow at Reynolds 1600. The numerical simulations, carried out with an accurate divergencefree scheme, outperform classical largeeddies formulations and provides a simple demonstration of the pertinence of the proposed largescale modeling.
skh]S. Kadri Harouna em]E. Mémin
Stochastic representation of the Reynolds transport theorem: revisiting largescale modeling
^{a}Laboratoire Mathématiques, Image et Applications (MIA), Université de La Rochelle, Avenue Michel Crépeau 17042 La Rochelle, France
^{b}INRIA Rennes, Campus Universitaire de Beaulieu, 35042 Rennes cedex, France
Keywords: Largescale fluid flow dynamics, stochastic transport, Subgrid model, turbulence,TaylorGreen flow
1 Introduction
The largescale 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 (Boussinesq77, ). 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 largescale dynamics in which the contribution of unresolved smallscale processes are explicitly taken into account. For the NavierStokes 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 partiallyknown 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.
Smallscale processes are responsible both for an energy dissipation but also for local backscattering of energy (Piomelli91, ). 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 (Leith90, ; Mason92, ; Schumann95, ). Recently those models have regained a great interest for the modeling of geophysical flows (Buizza99, ; Grooms13, ; Majda99, ; Majda03, ; Shutts05, ) in climate sciences (see also the thematic issue (Palmer08, ) or the review (Franzke15, )).
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 (Deardorff70, ; Lilly66, ; Meneveau00, ; Smagorinsky63, ; Sagaut05, ). This concept dates back to the work of Boussinesq (Boussinesq77, ) and Prandtl (Prandtl25, ). 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 timescale component and a fast oscillating velocity field representing the unresolved flow. Though there is, in general, no sharp timescale separation in turbulent flows, the resolved velocity can be interpreted as a temporally coarsegrained component whereas the timeuncorrelated component stands for the small timescale 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 (Ghosal95, ; Ghosal96, ). Besides, those equations introduce an effective advection related to the smallscale velocity inhomogeneity. This modified advection, empirically introduced in Langevin models of particle dispersion (Macinnes92, ), corresponds exactly to a phenomenon, termed turbophoresis, related to the migration of inertial particles in regions of lower turbulent diffusivity (Sehmel70, ).
The largescale representation of the NavierStokes equations on which we rely in this study are built from a stochastic version of the Reynolds transport theorem (Memin14, ). This modified Reynolds transport theorem, which constitutes here the cornerstone of our largescale dynamics representation, is presented in the following section. General invariance properties of the corresponding largescale 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 GreenTaylor vortex flow (TaylorGreen, ). We will show that all the proposed schemes outperform the usual dynamic Smagorinsky subgrid formulation (Germano91, ; Germano92, ; Lilly92, ; Smagorinsky63, ).
2 Stochastic modeling of fluid flow dynamics
Numerous methodological choices can be envisaged to devise stochastic representations of the NavierStokes 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 (BensoussanTemam73, ). Another choice, in the wake of Kraichnan’s work (Kraichnan59, ), consists in closing the largescale flow representation in the Fourier space by relying on a Langevin equation (Kraichnan70, ; Laval06, ; Leith71, ). 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 (Laval06, ; Schumann95, ). Lagrangian stochastic models based on Langevin equations have been also intensively used for turbulent dispersion (Sawford86, ) or in probability density function (PDF) modeling of turbulent flows (Haworth86, ; Pope00, ). Those Lagrangian models, which require to model the drift and diffusion functions, lead to very attractive particle based representations of complex flows (Pope94, ; Minier14, ). They are nevertheless not adapted to global largescale 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 (Memin14, ). 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,
(1) 
with the velocity components:
(2) 
In this decomposition the first righthand term is a smooth function of time associated to the largescale velocity component. The second term stands for the smallscale velocity field. It is a white noise velocity component defined from the (formal) timederivative of the random field: . This random field is a threedimensional 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 smallscale component is defined as a divergencefree random field; it is hence associated to a divergencefree diffusion tensor:
(3) 
Analogously to the standard deterministic case, the derivation procedure from the physical conservation laws of the NavierStokes equations is based primarily on the Reynolds transport theorem (RTT).
2.1 Stochastic 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 (2) with an incompressible smallscale velocity component (), this expression derived in (Memin14, ) (see also appendix Appendix A where a comprehensive derivation is provided for readers convenience), is given by:
(4) 
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 (Sagaut05, ) 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 smallscale component through the variance tensor, , defined from the diagonal of the smallscale 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 (Bardina80, ; Gent90, ; Lilly92, ; Smagorinsky63, ) or to the spectral vanishing viscosity (Karamanos00, ; Pasquetti06, ; Tadmor89, ). It is also akin to numerical regularization models considered in implicit models (Aspden08, ; Boris92, ; Dairay16, ; Lamballais11, ). 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 (Brooke92, ; Caporaloni75, ; Sehmel70, ). 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 (Macinnes92, ).
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 ():
(5) 
Together with the incompressibility of the random term, the incompressibility condition reads thus:
(6) 
In the case of an incompressible largescale flow component, , this reduces to:
(7) 
Note that for a divergencefree isotropic random field such as the Kraichnan model (Kraichnan68, ) the last condition is naturally satisfied, since this unresolved velocity component is associated with a constant variance tensor.
2.2 Mass conservation
Applying the RTT to the fluid density, , and canceling this expression for arbitrary volumes, we get the following mass conservation constraint:
(8)  
(9) 
For an incompressible fluid with constant density, together with a volumepreserving (isochoric) condition on the largescale velocity component, we retrieve the incompressibility conditions (7). It can be noted also that 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 smallscale component is exactly compensated by the energy dissipation associated to the diffusion term (Resseguier16a, ).
2.3 Linear 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 NavierStokes equations (Memin14, ):
(10a)  
(10b)  
(10c) 
In this expression is the dynamic viscosity, denotes the largescale (slow) pressure contribution and is a zeromean turbulent pressure related to the smallscale 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 largescale velocity.
It can be observed that the largescale 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 divergencefree isotropic random field (with a constant diagonal variance tensor), this system boils down to an intuitive constant eddy viscosity diffusivity model:
(11) 
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 smallscale 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 smallscale random velocity to live on the isodensity surfaces provides immediately a clear justification of the isopycnal diffusion employed in oceanic circulation models (Memin14, ). 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 datadriven strategies.
General invariance properties of the proposed largescale representation are listed in the appendix 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 quasinormal approximations (MoninYaglomB, ; Orszag70, ). 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 KolmogorovRichardson scaling and a dimensionless number relating the largescale kinetic energy to a characteristic value of the velocity variance at the resolved scale.
Interested readers may also consult (Resseguier16a, ; Resseguier16b, ; Resseguier16c, ) for the derivation of stochastic representation of several geophysical flows.
3 Numerical simulation and assessment
In this section we assessed, through numerical simulations, the performances of the proposed largescale 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 (Germano91, ; Lilly92, ). Those numerical experiments have been performed on the TaylorGreen flow (TaylorGreen, ).
3.1 TaylorGreen vortex flow simulation
TaylorGreen 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 smallscale 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 highorder 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 smallscale eddies due to vorticity increase and vortex stretching mechanism (Brachet83, ; Orszag74, ; TaylorGreen, ).
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 highorder numerical methods have been tested for the direct numerical simulation (DNS) of the TaylorGreen vortex flow, see (Chapelier13, ; Rees11, ) and references therein. In this study we used a discrete scheme built on a divergencefree wavelet basis (KadriIJCV13, ; Kadri13, ). This scheme presents several computational advantages. First of all, it guarantees a divergencefree solution in the physical domain with a good numerical complexity. Besides, as the spatial filters considered corresponds to a multiresolution 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 pseudospectral method.
3.2 Analysis 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
(12)
where is the rate of strain tensor: . The mean kinetic energy is linked to the dissipation rate by:
(13) 
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.3 Direct Numerical Simulation of the TaylorGreen flow
We evaluated first the ability of the divergencefree wavelet based method to reproduce faithfully the main features of the TaylorGreen 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 divergencefree wavelet has been run on a regular grid of points. For comparison purpose we performed a classical pseudospectral 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 (Chapelier13, ). 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 (Brachet83, ), numerous numerical experiences have been conducted to understand accurately the dynamic of the TaylorGreen vortex flow at this Reynolds number. As mentioned in (Chapelier13, ; Rees11, ), 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 divergencefree wavelet based simulation on a mesh grid. We also performed a Fourier pseudospectral 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 pseudospectral solution (with dealiasing) available from (Chapelier13, ). 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 divergencefree 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 divergencefree wavelet based method provides results of comparable quality to those obtained from pseudospectral 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.4 Variance 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 smallscale velocities to live on specific isosurface. Another route would consist in specifying the turbulent velocity components from smallscale measurements statistics. Despite all those directions are worth exploring as they open new strategies to design subgrid tensors, in this study we choose to focus on several simple models of tensor to explore the potential of the proposed modeling.
3.4.1 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 Appendix C.1.
3.4.2 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 multiscale projector. More precisely, since the velocity can be decomposed as:
(14) 
with the projector onto the scaling functions basis (Deriaz06, ; Kadri13, ), the resolved filtered velocity is defined by:
(15) 
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:
(16) 
The stochastic NavierStokes equations (10a) rewritten respectively for and reads:
(17) 
and
(18) 
Taking the difference of the momentum equations (17) and (18) for the two subsequent levels and provides the residual dynamics:
(19) 
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 lefthand side cancels after filtering. Then, instead of solving the righthand 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 quasinewton method combined with the BroydenFletcherGoldfarbShanno (BFGS) method Nocedal99 () 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 Appendix D.
3.4.3 Smagorinsky subgrid model
The third model evaluated corresponds to the classical Smagorinsky eddy viscosity formulation,
coupled with the Germano estimation procedure (Germano91, ; Lilly92, ) to fix dynamically the eddy viscosity constant from a least squares estimation (Lilly92, ) 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.5 Largescale simulation numerical results
The performances of the different subgrid 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 NavierStokes equation started from the largescale (filtered) initial conditions. It is thus not relevant to correlate the largescale 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 (DynS) 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 (OptSS) and the variance tensor built from the velocity temporal covariance (temporal). For both groups, we displayed on the righthand side of figure 4 and figure 5 the energy dissipation rate computed from the kinetic energy rate of change. On the lefthand 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 largescale 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 8 and 9. Wavelet power spectrum corresponds to an averaged version of the Fourier spectrum (Bruneau02, ) and exhibits the same slopes as the Fourier spectrum (Perrier95, ). 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 smallscale 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 (Germano91, ).
Regarding the different criteria explored, we see that the proposed largescale 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 twoscale 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 stateoftheart solutions for this flow (Chapelier13, ; Rees11, , and references therein). Though very simple, the temporal covariance model weighted by a unique constant fixed through the quadratic scaling rule presented section Appendix C.1 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. Largeeddies simulation models could likely be proposed in this spirit for nonhomogeneous 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.
4 Conclusion
In this paper we have described a decomposition of the NavierStokes 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 largescale representation paradigm, which can be interpreted as a large eddies simulation formalized through a timescale 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 largescale representation have been assessed on a TaylorGreen 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 smallscale observations. Another route of investigation consists in defining adapted basis for the smallscale random field from the fluctuations observed on two consecutive scales of the resolved tensor.
References
 [1] A. Aspden, N. Nikiforakis, S. Dalziel, and J.B. Bell. Analysis of implicit les methods. App. Math. Comp. Sci., 3(1):103–126, 2008.
 [2] J. Bardina, J.H. Ferziger, and W.C. Reynolds. Improved subgrid scale models for large eddy simulation. volume 80, page 1357.
 [3] A. Bensoussan and R. Temam. Equations stochastique du type NavierStokes. J. Funct. Anal., 13, 1973.
 [4] J.P. Boris, F.F. Grinstein, E.S. Oran, and R.L. Kolbe. New insights into largeeddy simulation. Fluid Dynamics Research, 10:199–228, 1992.
 [5] J. Boussinesq. Essai sur la théorie des eaux courantes. Mémoires présentés par divers savants à l’Académie des Sciences, 23 (1): 1–680, 1877.
 [6] M. Brachet, D. Meiron, S. Orszag, B. Nickel, R. Morf, and U.Frisch. Smallscale structure of the taylorgreen vortex. J. Fluid. Mech., 130(6):411–452, 1983.
 [7] M. Brooke, K. Kontomaris, T. Hanratty, and J. McLaughlin. Turbulent deposition and trapping of aerosols at a wall Phys. of Fluids, A 4, 825â834,1992.
 [8] ChH. Bruneau, P. Fischer, Z. Peter, and A. Yger Comparison of numerical methods for the computation of energy spectra in 2D turbulence. Part I: Direct methods Sampling Theory Sign. and Im. Proc., 1(1):0:50, 2002.
 [9] R. Buizza, M. Miller, and T.N. Palmer. Stochastic representation of model uncertainties in the ECMWF ensemble prediction system. Quarterly Journal Royal Meteorological Society, 125:2887–2908, 1999.
 [10] T. ChaconRebollo and R. Lewandowski. Mathematical and numerical foundations of turbulence models and applications. Springer NewYork, 2014.
 [11] J.B. Chapelier, M. De La Llave Plata, F. Renac, and E. Martin. Final abstract for ONERA Taylor Green DG participation. In 1st International Workshop on HighOrder CFD Methods at the 50th AIAA Aerospace Sciences Meeting, Nashville,TN, 2013.
 [12] M. Caporaloni, F. Tampieri, F. Trombetti, and O. Vittori Transfer of particles in non isotropic air turbulence. J. Atmos. Sci, 32:565â568,1975.
 [13] T. Dairay, E. Lamballais, S. Laizet, and J.C. Vassilicos. Numerical dissipation vs. subgridscale modeling for large eddy simulation. to be published, 2016.
 [14] J. Deardorff. A numerical study of threedimensional turbulent channel flow at large reynolds numbers. J. Fluid Mech., 1970.
 [15] E. Deriaz and V. Perrier. Divergencefree and curlfree wavelets in 2d and 3d, application to turbulent flows. J. of Turbulence, 7(3):1–37, 2006.
 [16] C. Franzke, T. O’Kane, J. Berner, P. Williams, and V. Lucarini. Stochastic climate theory and modeling. Wiley Interdisciplinary Reviews: Climate Change, 6(1):63–78, 2015.
 [17] P. Gent and J. McWilliams. Isopycnal mixing in ocean circulation models. J. Phys. Oceanogr., 20:150–155, 1990.
 [18] M. Germano. Turbulence : the filtering approach. J. Fluid Mech., 1992.
 [19] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamic subgridscale eddy viscosity model. Phys. of Fluids, 3:1760–1765, 1991.
 [20] S. Ghosal. An analysis of numerical errors in largeeddy simulations of turbulence. J. Comp. Phys., 125:187–206, 1996.
 [21] S. Ghosal and P. Moin. The basic equations for the large eddy simulation of turbulent flows in complex geometry. J. Comp. Phys., 1995.
 [22] A. Gronskis, D. Heitz, and E. Mémin. Inflow and initial conditions for direct numerical simulation based on adjoint data assimilation. J. Comp. Phys, 242(6):480–497, 2013.
 [23] I. Grooms and A. Majda. Efficient stochastic superparameterization for geophysical turbulence. PNAS, 110(12), 2013.
 [24] D. Haworth and S. Pope. A generalized Langevin model for turbulent flows. Phys. of Fluids, 29:387–405, 1986.
 [25] B. Dubrulle J.P Laval and J.C. McWilliams. Langevin models of turbulence: Renormalization group, distant interaction algorithms or rapid distortion theory? Pys. of Fluids, 15(5):1327–1339, 2006.
 [26] S. KadriHarouna, P. Dérian, P. Héas, and E. Mémin. Divergencefree wavelets and high order regularization. International Journal of Computer Vision, 103(1):80–99, 2013.
 [27] S. KadriHarouna and V. Perrier. Effective construction of divergencefree wavelets on the square. J. of Computational and Applied Math., 2013.
 [28] G. Karamanos and G. Karniadakis. A spectral vanishing viscosity method for largeeddy simulations. Journal of Computational Physics, 163(1):22–50, 2000.
 [29] R. Kraichnan. The structure of isotropic turbulence at very high reynolds numbers. J. of Fluids Mech., 5:477–543, 1959.
 [30] R. Kraichnan. Smallscale structure of a scalar field convected by turbulence. Phys. of Fluids, 11:945–963, 1968.
 [31] R. Kraichnan. Convergents to turbulence functions. J. of Fluid Mech., 41:189–217, 1970.
 [32] H. Kunita. Stochastic flows and stochastic differential equations. Cambridge University Press, 1990.
 [33] E. Lamballais, V. Fortunè, and S. Laizet. Straightforward highorder numerical dissipation via the viscous term for direct and large eddy simulation. J. Comp. Phys., 230:3270–3275, 2011.
 [34] W. Layton. Approximating the larger eddies in fluid motion v: Kinetic energy balance of scale similarity models. Math. and Comp. Modelling, 2000.
 [35] C. Leith. Atmospheric predictability and twodimensional turbulence. J. Atmos. Sci, 28, 1971.
 [36] C. Leith. Stochastic backscatter in a subgridscale model: plane shear mixing layer. Phys. of Fluids, 2(3):1521–1530, 1990.
 [37] D. Lilly. On the application of the eddy viscosity concept in the inertial subrange of turbulence. Technical Report 123, NCAR, 1966.
 [38] D. Lilly. A proposed modification of the Germano subgridscale closure. Phys. Fluids, 3:2746–2757, 1992.
 [39] J. Macinnes and F. Bracco. Stochastic particles dispersion modelling and the tracerparticle limit. Phys. Fluids, 4(12):2809–2824, 1992.
 [40] A. Majda, I. Timofeyev, and E. Vanden Eijnden. Models for stochastic climate prediction. PNAS, 1999.
 [41] A. Majda, I. Timofeyev, and E. Vanden Eijnden. A systematic strategies for stochastic mode reduction in climate. Journ. Atmos. Sci., 60:1705–1722, 2003.
 [42] P.J. Mason and D.J. Thomson. Stochastic backscatter in largeeddy simulations of boundary layers. J. of Fluid Mech., 242:51–78, 1992.
 [43] E. Mémin. Fluid flow dynamics under location uncertainty. Geophysical & Astrophysical Fluid Dynamics, 108(2):119–146, 2014.
 [44] C. Meneveau and J. Katz. Scaleinvariance and turbulence models for largeeddy simulation. Annu. Rev. Fluid. Mech, 32:1–32, 2000.
 [45] J.P. Minier S. Chibbaro and S. Pope. Guidelines for the formulation of Lagrangian stochastic models for particle simulations of singlephase and dispersed twophase turbulent flows. Phys. of Fluids, 26,1113303, 2014.
 [46] A.S. Monin and A.M. Yaglom. Statistical fluid mechanics, volume II. MIT Press, 1975.
 [47] J. Nocedal and S.J. Wright. Numerical optimization. Springer Series in Operations Research. SringerVerlag, NewYork, 1999.
 [48] S. Orszag. Analytical theories of turbulence. J. Fluid Mech., 41:363–386, 1970.
 [49] S. Orszag. Numerical simulation of the taylorgreen vortex. In SpringerVerlag, editor, International Symposium on Computing Methods in Applied Sciences and Engineering, volume 2, pages 50–64, 1974.
 [50] T. Palmer and P. Williams. Theme issue ’stochastic physics and climate modelling’. Phil. Trans. R. Soc., 366(1875), 2008.
 [51] R. Pasquetti. Spectral vanishing viscosity method for largeeddy simulation of turbulent flows. J. Sci. Comp., 27(13):365–375, 2006.
 [52] V Perrier, T. Philipovitch, and C. Basdevant,. Wavelet spectra compared to Fourier spectra J. Math. Phys., 36:1506â1519, 1995.
 [53] U. Piomelli, W. Cabot, P. Moin, and S. Lee. Subgridscale backscatter in turbulent and transitional flows. Phys. Fluids, 3(7):1766–1771, 1991.
 [54] S. Pope. Lagrangian PDF methods for turbulent flows. Annu. Rev. Fluid Mech., 26:2363, 1994.
 [55] S. Pope. Turbulent flows. Cambridge University Press, 2000.
 [56] L. Prandtl. Bericht uber untersuchungen zur ausgebildeten turbulenz. Angew. Math, Meth., 5:136–139, 1925.
 [57] G. Da Prato and J. Zabczyk. Stochastic equations in infinite dimensions. Cambridge University Press, 1992.
 [58] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part I Random transport and general models https://hal.inria.fr/hal01391420, 2016.
 [59] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part II Quasigeostrophy and efficient ensemble spreading https://hal.inria.fr/hal01391476, 2016.
 [60] V. Resseguier, E. Mémin, and B. Chapron. Geophysical flows under location uncertainty, Part III SQG and frontal dynamics under strong turbulence conditions https://hal.inria.fr/hal01391484, 2016.
 [61] P. Sagaut. Largeeddy simulation for incompressible flow  An introduction, third edition. SpringerVerlag, Scientic Computation series, 2005.
 [62] B. Sawford. Generalized random forcing in randomwalk models of turbulent dispersion models. Phys. of Fluids., 29:3582–3585, 1986.
 [63] U. Schumann. Stochastic backscatter of turbulence energy and scalar variance by random subgridscale fluxes. Proc. R. Soc. Lond. A, 451:293–318, 1995.
 [64] G. Sehmel. Particle deposition from turbulent air flow. J. Geophys. Res., 75:1766–1781, 1970.
 [65] G. Shutts. A kinetic energy backscatter algorithm for use in ensemble prediction systems. Quarterly Journal of the Royal Meteorological Society, 612:3079–3012, 2005.
 [66] J. Smagorinsky. General circulation experiments with the primitive equation: I. the basic experiment. Monthly Weather Review, 91:99–165, 1963.
 [67] E. Tadmor. Convergence of spectral methods for nonlinear conservation laws. SIAM J. Numer. Anal, 26(1):30–44, 1989.
 [68] G. Taylor and A. Green. Mechanisms of the production of small eddies from large ones. Proceedings of the Royal Society of London A, 1938.
 [69] W.M. Van Rees, D.I. Leonard, A.and Pullin, and P. Koumoutsakos. A comparison of vortex and pseudospectral methods for the simulation of periodic vortical flows at high reynolds numbers. J. Comp. Phys, 230:2794–2805, 2011.
Appendix Appendix A Stochastic 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:
(A.1) 
The random field involved in this equation is defined for all points of the fluid domain, through the kernel associated to the diffusion operator
(A.2) 
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:
(A.3) 
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 socalled quadratic variation process:
The notation stands for the quadratic crossvariation 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 Appendix B.
The drift term, , of Lagrangian expression (A.1), represents the ”smooth” – differentiable – part of the flow whereas the random part,
(A.4) 
figures the smallscale velocity component associated to a much thinner timescale. This component is modeled as a deltacorrelated random field in time as it represents a highly non regular process at the resolved time scale. In this work, we assume that this smallscale 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 (A.1) 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
(A.5) 
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 ItoWentzell formula [see theorem 3.3.1, 32]. This extended Ito formula incorporates in the same way as the classical Ito formula a quadratic variation terms related to process but also covariation terms between and the gradient of the random function . Its expression is in our case given by
(A.6) 
It can be immediately checked that for a deterministic function, the standard Ito formula is recovered since the covariations terms between and cancel.
It follows from A.6 that for a fixed grid point, function is solution of a stochastic differential equation of the form
(A.7) 
where the Brownian term must compensate the Brownian term of (A.6). The quadratic variation term involved in (A.6) is given through (appendix Appendix B) as
(A.8)  
(A.9) 
and the covariation term reads
(A.10) 
In these expressions (resp. ) designates the kernel associated to operator (resp. ). Now, identifying first the Brownian term and next the deterministic terms of equations (A.6) and (A.7), we infer
and finally get
(A.11)  
(A.12) 
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 (A.11) 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 righthand term of A.12 is selfadjoint) and if , where stands for the characteristic function, we get the sought form of this differential:
Appendix Appendix B Quadratic variation and covariation
We recall here the notion of quadratic variation and covariation, which play a central role in stochastic calculus. We will here restrict ourselves to standard Ito processes. Quadratic variation and covariation correspond respectively to the variance and covariance of the process increments along time. The quadratic covariation 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.
Appendix Appendix C General properties of the largescale stochastic NavierStokes model
First of all, it is important to emphasize that the distribution of the velocity anomaly, , with the Eulerian velocity field
(C.1) 
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 fourthorder correlation functions as considered in the Millionschikov hypothesis [46]. Nor is it based on quasinormal 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 NavierStokes equations. In the following sections we list several properties of this model. We begin first by a useful scaling relation of the variance tensor.
Appendix C.1 Variance tensor scaling
In the conditions of KolmogorovRichardson scaling, at scale , the velocity increments and the eddies turnover time scale as and respectively, with denoting a constant energy dissipation rate across the inertial scales range. The turnover 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:
(C.2) 
Within a given cell, , at scale , the energy of the smallscale random field is given as:
(C.3) 
which from (C.2) gives us:
(C.4) 
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.
Appendix C.2 Adimensionalization of largescale stochastic NavierStokes equations
Considering the scaled coordinates , , with velocity , pressure and a variance tensor , we have:
(C.5) 
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 largescale stochastic NavierStokes equations reads:
(C.6) 
where we introduced a dimensionless number
(C.7) 
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 (C.4) a characteristic value of the resolved variance tensor
(C.8) 
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, (C.6) tends to the original NavierStokes 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 largescale model proposed conserves the invariance properties of the original NavierStokes equations.
Appendix C.3 Invariance of largescale stochastic NavierStokes equations
Let us consider in the following the classical symmetry groups of the incompressible NavierStokes equations in a periodic domain and in the absence of external forcing.
Appendix C.3.1 Translationinvariance
Space translation invariance is achieved if the whole velocity field is still a solution of the NavierStokes equations. In our setting we have
Separating the Brownian component from the smooth term yields: , and hence