Development of stochastic models for turbulence

Subgrid scale interactions in forced homogeneous turbulence are decomposed into components correlated and uncorrelated with the resolved velocity field. The correlated part is well known to be equivalent to a spectral eddy viscosity. The mean energy transfer can be predicted even if the uncorrelated remainder is ignored, but statistics pertaining to the fluctuations of the nonlinear term require a model for the remainder, which is found to have a strongly non-Gaussian probability density. Some implications of including these non-Gaussian properties in stochastic models are described.


Introduction
Stochastic models have had an important role in turbulence theory in demonstrating the realizability of spectral closure theories (Kraichnan 1971, Kraichnan andLeith 1972). Recent theoretical developments of Langevin models for turbulence (Laval and Dubrulle 2006, Laval et al 2003 have been motivated in part by the possibility, suggested by Kraichnan (1976), of representing subgrid interactions by a stochastic model. The expectation is that models that include a definite representation of the subgrid scales will be more accurate than standard models that attempt to close subgrid interactions in terms of the resolved field alone. The present work is motivated by Kraichnan's (1985) proposal to construct stochastic models of subgrid interactions by imposing statistical constraints. Briefly, the subgrid interactions are replaced by a random force which is constrained to reproduce selected statistics of the actual interactions. It should be emphasized that this type of modelling does not aim at a practical model, since the constraints are evaluated from the exact solution; the goal instead is to understand the effect of replacing the exact interactions by a simpler stochastic surrogate.
This approach connects the formulation of a model to the properties of turbulence that the model is expected to reproduce; it suggests formulating a hierarchy of models that can have as many properties of Navier-Stokes turbulence as are considered relevant. For example, whereas a fairly simple spectral eddy viscosity model can predict the correct energy flux, a more complex model will be needed to predict energy transfer fluctuations and other aspects of 'intermittency' (Monin andYaglom 1975, Frisch 1995) in turbulence.
Whereas Kraichnan proposed to formulate the constraints using statistical symmetries, they are developed here directly from direct numerical simulation (DNS) data. The scope of this paper is limited to considering how Kraichnan's (1985) approach might be applied using numerical databases to construct the relevant constraints. The actual formulation and implementation of models is left for future research.

Spectral eddy viscosity models
Write the forced Navier-Stokes equations for homogeneous turbulence in the standard Fourier formu and P imn (k) = k m δ in + k n δ im − 2k −2 k i k m k n is the nonlinear coupling following elimination of the pressure. The random force f in (1) is concentrated at large scales and maintains the velocity field in a statistically steady state. Consider an arbitrary cutoff wavenumber k c and introduce the decomposition where so that N < only contains interactions among modes satisfying k = |k| k c , and all interactions with modes such that k k c occur in N > . Subgrid scale modelling replaces N > by some surrogateN so that the equatioṅ for the resolved velocity field u(k, t), which is defined on the smaller set of modes satisfying k k c , reproduces relevant statistical properties of the original velocity field satisfying (1): a typical minimum requirement is that the energy spectra computed from (1) and (5) coincide when k k c , but the possibility of reproducing higher order statistics can be considered as well. An interesting theoretical question is whether the stochastic process defined by (1), restricted to k k c , could be reproduced exactly (as a stochastic process, not necessarily for every realization) by a solution of (5) for a suitableN; this problem is addressed by Eyink (1996). The simplest model is the direct mode truncation N > ≡ 0. In forced isotropic turbulence, unless k c is very close to the Kolmogorov scale, the resolved field will simply grow in time under the forcing because energy cannot be removed sufficiently rapidly by viscosity alone. Equivalently, if the force is turned off, the energy in the large scales will be driven to equipartition instead of decaying. This trivial model must be rejected.
A general approach to stochastic subgrid modelling is suggested by Kraichnan (1985): the idea is to require that statistics ofN and joint statistics ofN and the resolved velocity field u agree with the corresponding statistics of N > . By requiring agreement of more statistics, a hierarchy of modeled velocity fields can be generated that share some, but not all, properties of Navier-Stokes turbulence. A major purpose of (Kraichnan 1985) is to describe how stochastic processes satisfying such (possibly nonlinear) constraints can be realized. We again stress that such models will have a theoretical rather than a practical value, since the constraints are formulated in terms of the known exact N > ; the goal is understanding of exactly what is lost in the replacement of the exact dynamics by a surrogate constructed by stochastic synthesis.

DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
A simple example is the constraint If it is satisfied, the steady inertial range energy balance of (1): is formally preserved by the model (5). The constraint (6) can be imposed directly onN following (Kraichnan 1985), but it also can be satisfied by a simple linear regression model: evaluate the correlation coefficient between u and N > : and setN This is the model considered by Domaradzki et al (1987) and more recently as the 'linear estimation model', a special case of a more general model, by Moser (1999, 2004). Other recent work introducing spectral eddy viscosity concepts includes (Bou-Zeid et al 2005; Laval and Dubrulle 2006). For additional applications of a decomposition of the unresolved nonlinearity into a components correlated and uncorrelated with the resolved velocity, see (Domaradzki and Loh 1999, Domaradzki and Saiki 1997. We consider some properties of this model. Data from 128 3 forced steady state DNS was used to evaluate the factor η(k|k c ) in (8). The corresponding eddy viscosity ν(k|k c ) defined by η(k|k c ) = ν(k|k c )k 2 is shown in figure 1 as a function of k for various values of the cutoff k c . The graphs are translated and rescaled linearly to a common abscissa. Consistently with many previous studies (Domaradzki et al 1987) and with notions of the locality of energy transfer, the viscosity has the 'plateau-cusp' structure predicted from closure analysis by Kraichnan (1976): it is constant for scales k k c well separated from the cutoff, and rises sharply near k c . Kraichnan's analysis applies to an infinite inertial range, giving the plateau limit ν(k|k c ) ∼ k −4/3 c for k k c and the precise behaviour near the cutoff ν(k|k c ) ≈ ν(k c |k c ) − (k c − k) 1/3 for k ↑ k c . Since the relatively low resolution DNS of the present and previous studies is very far from the conditions of this analysis, quantitative agreement with these predictions cannot be expected. Indeed, the predicted 'plateau' value is zero, and sometimes even slightly negative, except for the very low cutoff value k c = k max /8; again, similar results are found by Domaradzki et al (1987), Moser (1999, 2004), and in the context of scalar turbulence, by Kang and Meneveau (2005).
Explanations for the discrepancy have been proposed based on hypothetical departure of the spectral scaling from the Kolmogorov scaling assumed by Kraichnan (Chollet and Lesieur 1981); but since closure theories do not assume that the spectrum is Kolmogorov, or even that it obeys any power law scaling at all, we attempted to resolve the discrepancy by using a closure theory to compute the eddy viscosity using the measured energy spectrum. However, a preliminary calculation by the eddy-damped quasi-normal Markovianized (EDQNM) closure (Orszag 1973) did not result in significantly improved agreement. This issue remains an open problem for further investigation.  (8); the abscissa is the wavenumber rescaled and translated to be the same for each graph.

Models including a random force
It is natural to consider the modelling error defined by the 'remainder' With this definition, is merely a restatement of the exact dynamics; although this equation has the formal structure of a Langevin model, nothing suggests that R can reasonably be considered a 'random force'. Kraichnan's (1985) program can be continued by constructing modelsR for R in (11). Again, the simplest choice isR = 0, which is the model considered in the previous section: it is a spectral eddy viscosity model, comparable to the Chollet-Lesieur (1981) model, although the eddy viscosity is obtained from DNS data for N > rather than from a closure theory. It was noted earlier that (6) and (7) suggest that this model can preserve the energy balance of (1) and consequently preserve the energy spectrum for k k c ; this expectation is borne out by numerical tests (Langford and Moser 2004).
We pass to more complex models by asking how closeN = ηu is to N > . This question was addressed by Langford and Moser (1999) by evaluating the ratio R(k, t)R(−k, t) / N > (k, t)N > (−k, t) in (11): a small value indicates that the approximation is good, whereas a value near one indicates that the model estimate is poor. The corresponding data for our 6 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 2. Spectrum of the remainder term defined by (10) for cutoff wavenumber k c = k max /4 (curve with symbols). The straight line is the power law k 4 for comparison. The ordinate is the spectrum defined by approximate integration of the correlation R p (k)R p (−k) over a spherical shell of radius k = |k| and the abscissa is the wavenumber k.
computations (not shown) were in close agreement with those of (Langford and Moser 1999): the variance of N > nearly equals the variance of R = N > − ηu, essentially over the entire range of resolved wavenumbers. The result is obviously expected for k k c , when η ≈ 0, but even when k ≈ k c , the ratio was not close to zero.
These calculations demonstrate that the fluctuations of N > are not well modeled byN = ηu. A possible conclusion is that the model withR = 0 is simply inaccurate, and that accuracy requires some nontrivial model for the remainder R. But another view is that if the mean energy transfer can be predicted by assuming thatR = 0, it must be the case that most of the fluctuations of N > do not participate in mean energy transfer. This view is supported by numerical experiments which demonstrate that net energy transfer is the result of cancellations between large instantaneous fluctuations in forward and backward transfer (Gotoh and Watanabe 2005).
Nevertheless, if we want to model quantities related to fluctuations of N > , it is necessary to consider the statistical properties of R more closely. The most basic statistical property of R is its spectrum |k|=k R i (k, t)R i (−k, t) .
(12) Figure 2 shows a typical spectrum evaluated from the DNS data. Over a nontrivial range of scales, this spectrum has nearly a k 4 slope, in apparent agreement with closure analysis suggesting that R is related to the energy 'backscatter' from small to large scales (Lesieur 1987). The result for 7 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 3. The remainder spectrum of figure 2 for a range of cutoff wavenumbers k c : squares, k c = k max /8; triangles, k c = k max /6; deltas, k c = k max /4; diamonds, k c = k max /3; circles, k c = k max /2. The ordinate is as defined in figure 2, but the abscissa is the wavenumber axis translated so that all graphs end at the same point. different cutoffs appears in figure 3; the graphs are shifted so that the cutoff appears at the same abscissa value. As expected, the force amplitude increases as the cutoff k c approaches zero. The apparent agreement with closure must be qualified because the k 4 backscatter spectrum should apply to the largest scales of motion, not to intermediate scales. But in this case, it is important to consider the role of finite resolution, and in particular the inability of DNS forced at nearly the largest scales of motion, to treat asymptotically large scales.
We would like to briefly outline some possible modelling strategies for the remainder term. Replacing (11) bẏ defines a stochastic model with a random force of the type considered for example in (Chasnov 1991, Laval and Dubrulle 2006. The simplest model is a GaussianR chosen to have the same spectrum and/or time correlation structure as the actual remainder term R. It should be noted that time correlation issues cannot be avoided if we wish to propose a model for R: they were evaded in the eddy viscosity model by postulating the proportionality ofN and u. In this case, (13) is a system of nonlinear coupled Langevin equations for the Fourier amplitudes u i (k, t) from which a Fokker-Planck equation for the joint probability density function (pdf ) of these amplitudes can be derived. We next consider how we might formulate more refined statistical models based on the pdf of R. A first question is whether this pdf is Gaussian. Ideally, we would like to compute it at all wavenumbers; however, to obtain a statistically adequate sample size, the pdf of R i (x|k c ), obtained by transforming R i (k|k c ) to physical space, was evaluated instead. The individual 8 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Figure 4. Semilogarithmic plot of the pdf of the remainder in (10) after inverse Fourier transform to physical space; only the x-component is shown due to isotropy. The pdf of a Gaussian with the same mean and variance is shown (symbols) for comparison. The abscissa is the value of the fluctuation of the remainder; the ordinate is the relative frequency.
components R i (x|k c ) were all found to have a flatness R 4 i / R 2 i 2 ≈ 5, well in excess of the Gaussian value of 3. The measured pdf of the component R x is compared to a Gaussian pdf with the same mean and variance in figure 4: although the resolution inevitably degrades in the tails of the pdf, the characteristic 'heavy tails' (that is, higher than Gaussian probability of extreme events) are evident in the semi-logarithmic plot.
More precise data appear in table 1. The variance and flatness of the low-and high-pass filtered fields obtained by transforming H(k * − k)R i (k|k c ) and H(k − k * )R i (k|k c ) (as usual, H(x) = 0 for x < 0 otherwise H(x) = 1) to physical space, are shown for various choices of k * < k c . The variance values provide a check on the data processing because the sum of the variances of the low-and high-pass filtered fields must be independent of k * . Although the lowpass filtered data has much larger flatness than the high-pass filtered data when k * = 0.5k c and 0.75k c , the extremely low amplitude of the filtered field in these cases makes these statistics unreliable. What is perhaps surprising is that when k * = 0.925k c is extremely close to the cutoff and the variances in the high-and low-pass filtered fields are almost equal, the high-pass data has a flatness not far from the Gaussian value. Nevertheless, the large flatness in the low-pass filtered data, which can now be considered statistically significant, shows that the remainder field is robustly non-Gaussian.
We consider how these non-Gaussian properties could be incorporated in a modelR for R. One possibility is to defineR by a polynomial mapping of a correlated Gaussian field. A stochastic model of this type is considered by Luczka et al (1995). In the limit that the correlation time of the mapped Gaussian field is very small, the result is a general Markov model in which the pdf of 9 DEUTSCHE PHYSIKALISCHE GESELLSCHAFT Variance k k * Variance k k * Flatness k k * Flatness k k * k * = k c /2 4.60 × 10 −7 1.147 × 10 −4 25 5.2 k * = 3k c /4 5.58 × 10 −6 1.096 × 10 −4 16 4.6 k * = 0.925k c 4.05 × 10 −5 7.47 × 10 −5 7.6 3 .5 the velocity field satisfies an evolution equation which is nonlocal in k, rather than the Fokker-Planck partial differential equation that would result from GaussianR (Kubo 1985). While the nonlinear term N < can induce non-Gaussian properties of the modelled velocity field, these properties originate in the large scales; a model with non-Gaussian random forcing can include non-Gaussian small-scale effects as well.

More general subgrid models
We briefly note two alternatives to the Markovian linear regression expressed by (8) on which this discussion is based. The more general linear non-Markovian regression hypothesis leads to an integral equation for the damping, This hypothesis allows more realistic modelling of time-correlations and permits stochastic constraints based on two-time properties. Kraichnan (1985) considers this type of modelling exclusively. The resulting stochastic models are close in structure to the original generalized Langevin model for the direct interaction approximation (Kraichnan and Leith 1972). Another direction is to abandon linear regression altogether by decomposing N > into a part statistically dependent on u(k, t) and a remainder: where L > i (k|k c , t)|u(k, t) denotes the expectation of L > i (k|k c , t) conditioned on u(k, t). This decomposition is much more general than (8) because the conditional expectation need not be linear in u. The optimal LES of Langford and Moser (1999) is obtained by settingR = 0: single-point single-time moments of u would be predicted exactly by this model regardless of R (Pope 2000), although a model for the remainder term can be relevant to the two-point and two-time statistics of interest in the theory of homogeneous turbulence. Computations in which the conditional expectation is systematically approximated by a quadratic model show rather small departures from the linear model (Langford and Moser 1999); this observation implies that quadratically nonlinear modelling does not significantly alter the fluctuating N > field.

Conclusions
We have considered stochastic models of subgrid scale interactions in turbulence based directly on DNS data with a view toward applying the idea of realization under stochastic constraints introduced by Kraichnan (1985).
The simple linear constraint (6) leads to a spectral eddy viscosity model that can reproduce the mean energy flux. In agreement with previous studies based on the same idea of decomposing the unresolved nonlinear term into parts correlated and uncorrelated with the resolved velocity field (Chasnov 1991, Domaradzki and Loh 1999, Domaradzki and Saiki 1997 the eddy viscosity exhibits the 'plateau-cusp' structure shown in figure 1. As in numerous recent developments of Langevin models for turbulence , Laval and Dubrulle 2006, Laval et al 2003 we considered the possibility of stochastic models based on a model for the formal 'remainder' term R in (10). The energy spectrum of the remainder term shown in figures 2 and 3 provides the simplest statistical characterization of the fluctuation properties of the unresolved nonlinear term N > . As one example of turbulence properties related to fluctuations of N > , we mention the fluctuations in energy transfer, which have been measured by Gotoh and Watanabe (2005).
More refined properties of the fluctuations of the nonlinear term are provided by the pdf in figure 4 and the flatness values of table 1. Both point to significant departure from Gaussianity in the remainder term. These properties may be relevant to at least some kinds of 'intermittency' in turbulence. Development of stochastic models that can reproduce some of these non-Gaussian properties is an important topic for future research.