Statistical description of turbulent dispersion

• A submitted manuscript is the author's version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


I. INTRODUCTION
The statistical description of turbulent flow is known as one of the unsolved problems of physics. Only partial results exist, notably the theory of small viscous isotropic scales of turbulence of Kolmogorov [1,2]. A general theory of the statistical process of turbulence has yet to be developed [1].
One of the unsolved issues in statistical turbulence has been the description of dispersion. The conventionally adopted approach is to describe turbulent dispersion by a diffusion model [1]. But the model rests on empiricism and lacks a thorough mathematical basis [1,3,4]. Rigorous derivation of the diffusion equation is possible only for the case of homogeneous stationary turbulence [1,3,4]. Real turbulence is not of that kind. Applying the diffusion model in practical circumstances resorts to postulating diffusion coefficients on the basis of dimensional reasoning and fitting unknown constants to experimental results.
A more fundamental approach is that which starts from a Langevin equation for the velocity of a marked fluid particle or passive fluid admixture. The approach complies with the asymptotic structure of large-Reynolds-number turbulence according to Kolmogorov [1]. Many aspects of this theory have, in the meantime, been confirmed by experiment [1,5], theory [2], and direct numerical simulations (DNS) of the Navier-Stokes equations [6]. This includes the Lagrangian version of Kolmogorov theory, which is relevant for dispersion modeling. Studies published in the past decade confirm the universal statistical properties of Lagrangian small-scale turbulence [7][8][9]. One of the features is that for large Reynolds number Re = u * Lν −1 , where u * is the typical magnitude of the fluctuating fluid velocity, L is the length scale of the configuration of the fluid flow, and ν is the kinematic viscosity, accelerations of particles become more and more δ correlated compared to those of velocity [1]. It forms the basis for describing the statistics of fluid particle velocity by a Markov process. Required values of the Reynolds number Re for Kolmogorov theory to be applicable are about 10 3 and larger. This covers most practical applications of flows in nature and technical devices. An approach which starts from a Langevin model is, thus, from a physical and practical point of view, a most promising one.
Implementing the Lagrangian version of Kolmogorov's inertial subrange limit leads to a Langevin equation in which the white-noise term is of specified isotropic form. Specification of the Langevin equation can then be completed on determining the damping function. This has been successfully accomplished for homogeneous isotropic decaying turbulence [6,[10][11][12][13], the type of turbulence which occurs behind a grid of bars through which fluid with constant speed flows. But specification of the damping function for general inhomogeneous turbulent flow has been an issue which has yet to be resolved. Several of the proposals which have been made resort to empirical extensions to the exact descriptions for homogeneous isotropic turbulence [6].
An important step forward in the determination of the damping function was made by the introduction of the wellmixed principle [4]. The Fokker-Planck equation associated with the Langevin equation allows a Eulerian interpretation. The distribution of fixed-point fluid velocity should satisfy this equation. The Eulerian velocity distribution can be assumed to be known and the requirement that it should satisfy the Eulerian version of the Fokker-Planck equation constraints the form of the damping function. The requirement was first introduced in generalized form by Thomson [4]. It has since then been applied by various scientists in attempts to define the damping function for various sorts of turbulence [4,11,[13][14][15]. The general outcome of these efforts was that application of the well-mixed principle leads to complete specification of the damping function for homogeneous isotropic decaying turbulence only. Although the principle severely constraints the damping function for general anisotropic inhomogeneous turbulence, a general solution could not be derived. It is known as the nonuniqueness problem [11,[13][14][15].
A way out of the nonuniqueness problem was recently described in two subsequent papers: one focused on wallinduced turbulence [16] and the other on Gaussian anisotropic inhomogeneous turbulence [17]. The idea put forward was to exploit the largeness of the universal Kolmogorov constant C 0 present in the white-noise term of the Langevin equation. For values of the Reynolds number large enough for Kolmogorov theory to be applicable, the value of C 0 is about 6 [6,9,14,16,18]. For large values of C 0 , correlations of fluid particle velocity decrease in a short time. During 066309-1 1539-3755/2012/86(6)/066309 (14) ©2012 American Physical Society this short time, the change of energy because of dissipation will have hardly changed. The underlying mechanics of fluid particle motion will be close to that of a Hamilton process. Analogous to the methods applied in statistical mechanics of molecular motion [19][20][21], one can impose Onsager symmetry. It leads to complete specification of the first term of an expansion in powers of C 0 −1 of the damping function. Terms next to leading order can be specified as well. Moreover, expansion in powers of C 0 −1 also allows a derivation of the diffusion equation with Eulerian-based specification of the diffusion coefficient for Gaussian inhomogeneous anisotropic turbulence. The expressions are valid up to truncation of terms of O(C 0 −2 ). Objective of the present analysis is to deepen and generalize the mathematical and physical basis of the derived Langevin and diffusion equations. To do so, we start from a general Langevin model for fluid particle velocity which complies with the outcome of Kolmogorov's theory for large-Reynoldsnumber turbulence. From here onward we develop solutions solely by expanding in powers of C 0 −1 . We do not adopt any other a priori assumptions such as Gaussianity of the Eulerian flow field as done previously [16,17]. We determine the leading-order terms by employing the results of classical Hamiltonian-based statistical mechanics. We present nextto-leading-order descriptions which exhibit the features of non-Hamiltonian behavior related to change of energy on the time scale of velocity correlation. We assess the accuracy of the descriptions by comparison with results from experiments, theory, and DNS of decaying grid turbulence, pipe flow, and channel flow. General applicability of the presented descriptions is analyzed and discussed.

II. FORMULATION AND METHOD OF SOLUTION
The correlation time of the acceleration of a fluid particle can be represented by the Kolmogorov time scale. Compared to the correlation time of fluid particle velocity it decreases as Re −1/2 with increasing Reynolds number [1], a feature confirmed by experiment [5]. It gives rise to the idea to describe the statistical process of fluid particle velocity variation by a Langevin equation with the white-noise term described according to the inertial subrange of Lagrangian Kolmogorov theory [4,11,15]: where t is the time and v i = v i (t) is the statistical representation of fluctuating fluid particle velocity relative to Eulerian mean velocity u 0 i (x) evaluated at particle position x = x(t). Velocity is related to position by In Eq. (1), a i = a i (v ,x) is the unknown damping function, C 0 ≈ 6 the universal Kolmogorov constant, ε = ε(x) the energy dissipation rate averaged at fixed position, and w i (t) the Gaussian white noise of unit intensity. Considered is stationary turbulence of an incompressible fluid. Fixed-point statistical averages of Eulerian flow variables such as mean flow u 0 (x) and energy dissipation rate ε(x) are constant with respect to time and are assumed to be known. Statistical averages connected to marked fluid particles or passive admixture moving through the fluid are to be determined by solution of Eqs. (1) and (2). A consequence of the Markov model is that the intermittent behavior typical for small scales [8] has been disregarded. Such behavior is most apparent in the statistical averages of time and spatial derivatives of velocity but less pronounced in velocities itself and even small in displacements. Fractal models capturing intermittency in accordance with refined Kolmogorov hypothesis reveal small effects of such intermittency in the statistical distributions of particle dispersion [22]. Equation (1) therefore constitutes a solid basis for the statistical description of turbulent dispersion. At issue is the description of the damping function a (v ,x) as a function of particle velocity and position. It is this function which largely determines the statistics of dispersion. The aim and challenge are to determine the damping function and with it the Langevin equation for general turbulent flow. The subsequent aim is to derive the diffusion equation and its diffusion tensor D ij (x) from the Fokker-Planck equation associated with the Langevin equation.
To determine the Langevin and diffusion equations, we introduce a solution procedure in which we exploit the smallness of C 0 −1 . We introduce a perturbation scheme in which variables are described by expansions in powers of C 0 −1 . Equating terms of equal powers of C 0 −1 and solving leads to the desired descriptions for successive terms of a i and D ij (x). All this is rather standard from the viewpoint of regular perturbation techniques [23]. What differs is that, in the present case, the small parameter C 0 −1 is not a dimensionless number in the usual sense, viz. a dimensionless combination of physical quantities which can be varied in magnitude such that it can be made as small as one likes. Unfortunately, such a dimensionless combination cannot be identified in the present problem of the Langevin and Fokker-Planck equation for general inhomogeneous anisotropic turbulence: All terms in the Langevin and Fokker-Planck equation scale by the same combination of quantities representing velocity and length [17]. A way out is to exploit the presence of C 0 −1 to identify and order the magnitude of terms. Although C 0 is a constant, it can be treated as an autonomous scaling parameter because it appears as a separate entity in the statistical model. The Kolmogorov constant enters into the statistical model via the Lagrangian representation of the inertial subrange limit of small-scale turbulence. The statistics of the small scales are known to decouple from those of the large energetic scales in turbulent flow at large Reynolds number. Statistical quantities which appear as parameters in the Langevin and Fokker-Planck equations and which are governed by the large scales are independent from C 0 . It enables us to specify the dependency on C 0 of all terms in the Langevin and Fokker-Planck equations and we can use this information to set up a consistent and well-defined perturbation scheme based on expansion in powers of C 0 −1 . From a mathematical point of view any value of C 0 −1 is allowed. We can make C 0 −1 as small as we like to arrive at a situation where the first term of the expansion dominates the solution. But in turbulence C 0 −1 is limited to a value of about 1 6 . The question is whether such a value is small enough to make the expansion work. Do subsequent terms decrease 066309-2 sufficiently in magnitude and does the behavior revealed by the solutions comply with that observed in real turbulent flow? To deal with these issues we derive descriptions not only for terms of leading order in C 0 −1 (Sec. III) but also for terms of next to leading order in C 0 −1 (Secs. IV and V). We shall analyze the derived expressions in detail and compare them with results from theory, experiment, and DNS for representative cases of turbulent flow (Sec. VI and accompanying appendices). General conclusions regarding accuracy and applicability are made.

III. LEADING-ORDER DESCRIPTION
Attention is focused on passively marked fluid particles which are at time t = t 0 at position x = x 0 . The Kolmogorov constant is considered to be a scaling parameter whose inverse value is small. A meaningful description for the statistical velocity of the particles follows from Eq. (1) only if all terms are retained when collecting terms of leading order in C 0 −1 . That is, all three terms in the Langevin equation (1) should scale in the same manner with respect to C 0 . White noise of unit intensity scales as t −1/2 . Equating the orders of magnitude of the left-hand side of Eq. (1) and the white-noise term on the right-hand side then yields for time the scale C −1 0 u 2 * ε −1 , where u * is the typical magnitude of fluctuating velocity, i.e., standard deviation. The combination C −1 0 u 2 * ε −1 is known as the Lagrangian time scale [24]. The damping term in Eq. (1) equals in order of magnitude the two other terms if we scale a i by C 0 εu −1 * . In accordance with the above identified scalings with respect to C 0 we introduce, for t and a i , Langevin equation (1) then becomes For the position of the particle we can write where x * i is given by ] which do not depend on C 0 . Solutions of these equations will not depend on C 0 either. We can conclude that the process of statistical velocity takes place on a time scale of t * which does not depend on C 0 . In terms of the original time t this corresponds to a process that takes place on a time that is proportional to C −1 0 . Times by which particle velocities decorrelate scale in the same manner.
When following with time a fluid particle according to the descriptions of Eqs. (5) and (7), the coefficients a * i and ε will change in magnitude because of inhomogenity. But when C 0 1 the displacement of a particle during the time that there exists a correlation between particle velocities will be small to the extent that a * i and ε will have hardly changed in value compared to their value at x = x o . We can assume a * i and ε to be constant with respect to x * as long as we limit ourselves to the leading-order representation with respect to C 0 . This approximation also holds in the direction of the mean flow. The mean flow is generally large in turbulent flow, leading to large displacements per unit time. But changes in values of Eulerian-based statistical averages remain small. In a frame which moves with the mean flow, changes in the statistical averages generally take place over the external length scale L; L is equivalent to u 3 * ε −1 [1]. Displacement because of fluctuating velocities, however, is typically C −1 0 L during the time that there is correlation between velocities. It is small compared to the scale of inhomogenity when C 0 1. For more explicit analytical presentations of this approximation scheme, we refer to the cases of decaying grid turbulence [17] and Appendix A and the log-layer of wall-induced turbulence [16]. Now given that ε and a * i are constant with respect to x and equal to their values at x = x 0 , i.e., ε = ε 0 and a * i = a * i0 , the next question is how does a * i0 depend on the statistical velocity v i . For that purpose, we resort to the results of classical Hamiltonian-based statistical mechanics. Fluid flow is dissipative non-Hamiltonian. Energy is continuously pumped into the flow, which is dissipated through viscosity, From this relation it follows that the time scale of change of energy is u 2 * ε −1 0 while velocity decorrelation takes place on the time scale C −1 0 u 2 * ε −1 0 . The time scale over which the energy change takes place is, thus, a factor C 0 longer than the typical time of velocity fluctuations or the time over which particle velocities decorrelate. During the short time of velocity correlation, the mechanics of the velocity of the fluid particle are those of an almost Hamiltonian process, Several accounts exist on Hamiltonian treatments of inviscid flow [25]. But we do not have to resort to these: The knowledge that the flow can be assumed Hamiltonian is sufficient to arrive at statistical descriptions of the velocity process in the leading-order representation with respect to C 0 −1 . We thus can circumvent the complexities of the connections between descriptions at the micro and macro levels [26].
To establish the connection with classical statistical mechanics, we have chosen to consider the random motion of a fluid particle which can be thought of being in a box of linear dimensions l 1 ∼ LC −γ 0 , where L is external length scale of the flow, e.g., the radius of a pipe in which turbulent flow takes place, height of a channel, thickness of the atmospheric boundary layer, and so on. The box moves with the mean flow and its size is small to the extent that inside the box ε and a * i can be taken constant with respect to x and equal to their values at x = x 0 . As C 0 1 this is achieved by taking γ > 0. At the same time, the box is taken large enough such that the fluctuating velocity has established a complete statistical process before the randomly moving particle has left the box, i.e., one can construct a complete power-density spectrum from the time record established inside the box and the velocity has become uncorrelated from its initial value when leaving the 066309-3 box. The time that the particle is in the box is l 1 /u * . The characteristic time of the velocity process scales as C 0 −1 L/u * so decorrelation takes place inside the box as long as γ < 1. The physical laws which govern the motion of a particle in the box can be treated as Hamiltonian. In actual cases of turbulence there will be homogeneous directions in space and/or time, along which the statistical averages equal those at x = x 0 . An infinite number of boxes or realizations of the same statistical process thus can be identified in the configurations of turbulent fluid flow.
In statistical mechanics the distribution of velocity is usually taken Gaussian. Also in turbulent fluid flow this is an obvious choice. In many cases of turbulent fluid flow the fluctuating velocity is a small perturbation on a large mean flow. Under these circumstances one can adopt Einstein fluctuation theory [21] to show that the distribution of fluctuation velocity should approach that of a Gaussian distribution. In situations where the mean flow is not large, a Gaussian distribution will prevail when the Hamiltonian is quadratic in velocity [20,21]. The implication is that, in the leading-order formulation with respect to C 0 , the velocity distribution can be expected to be Gaussian. This conclusion is supported by a lot of experimental evidence [5,[27][28][29][30][31] and recent results of DNS [32]. They all point to Gaussian velocity distributions to a degree where skewness and kurtosis (kurtosis is flatness minus 3) are limited to values of about 0.3 and less. This would well comply with a perturbation expansion of the probability distribution in which the leading term is Gaussian and the next-to-leading-order where repeated indices imply summation. The fluctuation-dissipation theorem yields a connection between the damping term and the fluctuation term in Eq. (5) such that where the overbar denotes Lagrangian averaging, viz. the average value of many realizations at some fixed moment in time assuming that the particles are passively marked, i.e., the velocities of the particles have randomly chosen values at t = t 0 and x = x 0 in accordance with the statistical distribution of Eulerian-based velocities at the fixed point x 0 . As already mentioned, for we have a quasihomogeneous situation where fixed-point Eulerian-based statistical averages do not vary with space. Under homogeneous circumstances Eulerian and Lagrangian velocity distributions are the same [33]. That is, as long as t − t 0 = O(C 0 −1 ) one can take the statistical distribution for particle velocity equal to that of the fixed point Eulerian velocity at x = x 0 . Both are Gaussian.

Equality of Eulerian and Lagrangian velocity statistics implies
as long as t − t 0 = O(C 0 −1 ); σ ij 0 is the value of Eulerian covariance or Reynolds stress at x = x 0 . Hamiltonian dynamics are invariant to time reversal so Onsager symmetry holds [20,21], The fluctuation-dissipation relationship, the Euler-Lagrangian connection, and Onsager symmetry make that where λ ij 0 = σ −1 ij 0 , so the damping function in the Langevin equation becomes where we dropped the subscript 0 as the equation can be applied at any position, i.e., the values of λ ij = λ ij (x) and ε = ε(x) are adapted to their local values when moving through the inhomogeneous field while following a marked particle. The instantaneous adaptation is allowed as long as the damping is large, i.e., when C 0 1, so decorrelation of velocities occurs in locally homogeneous areas.
Apart from analogies, there are some striking differences of the above results with the statistical description of a particle in a bath of molecules as known from statistical mechanics. In statistical mechanics the damping function generally follows from a connection to macroscopic damping, such as Einstein's application of the Stokes' formula valid for damping of a large particle in a continuum of fluid [21]. The unknown coefficient of the white-noise term is subsequently determined using the fluctuation-dissipation theorem [34]. For a fluid particle in turbulent flow the opposite route is followed. The white-noise term is fully determined by the connection to the inertial subrange representation of Kolmogorov's universal similarity theory. The damping function is subsequently determined by applying the fluctuation-dissipation theorem supplemented with Onsager symmetry. [Note that the fluctuation-dissipation theorem in itself is not sufficient to specify α ij ; Onsager symmetry is also necessary to arrive at solutions (14) and (15).] The resulting expression for damping also allows a fluid mechanical interpretation. Therefore, it is noted that ε scales as u 3 * /L. One can then show that the damping term corresponds to the inviscid fluid force which acts on a "lump of fluid" or "agglomeration of fluid particles" of linear dimensions ∼LC 0 −1 . As C 0 1, this is small compared to the previously discussed box of dimensions LC −γ 0 , 0 < γ < 1. Explicit and general expressions for the damping force as solutions from the conservation laws of inviscid fluid flow have not been derived so far. What comes closest are the expressions for velocity-dependent lift force on a body in rotational inviscid flow [35].
Equations (1), (2), and (15) enable the simulation of particle tracks and the determination of their statistics. Analysis of particle dispersion is also possible from the diffusion equation, which is the mathematical approximation of Eqs. (1), (2), and (15) valid for long times after marking. The diffusion equation is most conveniently derived from the Fokker-Planck equation for joint probability p of velocity and position of fluid particle associated with Eqs. (1), (2), and (15), where t is time in a coordinate system which moves with the mean velocity, 066309-4 Solutions from Eq. (16) can be constructed through expansion in powers C 0 −1 ; see Refs. [16,17]. The method is analogous to the solution of Kramer's equation [34] and can be summarized as follows. Describe p by Substituting the above expansion in Eq. (16) and equating terms of equal power in C 0 −1 yields for p 0 and p 1 the equations Equation (19) has the solution where where λ is the determinant of λ ij and where G 0 = G 0 (x,t ) = ∞ −∞ pdv is probability density of particle position or admixture concentration. To determine G 0 we need to consider the descriptions of higher-order terms. Integrate Eq. (20) with respect to v over the entire v domain. The right-hand side must be zero because probabilities are assumed to go exponentially fast to zero as |v | → ∞. Substituting in the left-hand side solution (21) and noting that p G is zero mean, we find that (∂/∂t )G 0 = 0. In other words, G 0 must be constant on the time scale of t . Characteristic for the diffusion limit is that probability distributions of particle position vary over a much longer time scale. To model this we introduce The equation for p 2 then follows from Eq. (16) as The equation which determines G 0 = G 0 (x,t ) viz. the diffusion equation can now be derived from Eqs. (20) and (24) An expression for the right-hand side can be derived from Eq. (20). Multiply Eq. (20) with v k and integrate with respect to v over its entire domain. Applying partial integration to the term on the right-hand side, we obtain Multiplying with σ nk , which is the Eulerian covariance or Reynolds stress [cf. Eq. (12) σ nk = λ −1 nk ], and substituting solution (21) yields Substituting this result into Eq. (25), one obtains Returning to the time of the original nonmoving coordinate system [cf. Eqs. (17) and (23)], this becomes Equation (29) is the diffusion equation in the leading-order representation with respect to C 0 −1 . When σ nk is constant with respect to x we can combine the coefficients in the diffusion term to a single diffusion constant: 2C 0 −1 ε −1 σ in σ nk . For nonconstant σ nk a single diffusion constant occurs on adding the drift term (∂/∂x k )σ ik to damping function (15) and repeating the above solution procedure. The analysis of Sec. V will include the next-to-leading-order terms in the damping function and will yield a single diffusion coefficient.

IV. HIGHER-ORDER FORMULATION OF THE LANGEVIN EQUATION
So far attention has been focused on the leading-order term in the expansion with respect to C 0 −1 . The resulting descriptions for dispersion, i.e., Eqs. (1), (2), (15), and (29), will involve a truncation error of O(C 0 −1 ). Such an error will become smaller the larger C 0 is. But in turbulence the value of C 0 is limited to about 6. This corresponds to C 0 −1 = 0.17 and implies that the truncation error can become quite large. Deriving expressions for higher-order terms is, thus, wanted. For that purpose use is made of the well-mixed principle of Thomson [4]. The principle follows from the Eulerian interpretation of the Fokker-Planck equation for joint probability p of particle velocity v i and position x i associated with Eqs. (1) and (2), Given some initial distribution, particles will, in the course of time, mix up with the fluid and attain the distribution of fluid velocity. This equilibrium distribution satisfies the Eulerian interpretation of Eq. (30), which is given by where a i = a i (u ,x) is to be determined. The Eulerian velocity distribution can be taken Gaussian to leading order, where is the correction on Gaussian behavior. Values of zero, first-, and second-order moments are fully captured by the Gaussian part of the description, values of cumulants higher than second order are determined by f c (u ). Substituting the description for p E and Eq. (32) into Eq. (31), one obtains for a i the equation where we dropped terms of relative magnitude O(C 0 −1 ) in the contributions due to non-Gaussianity, i.e., the second term on the right-hand side of Eq. (35). Equation (35) is exact, i.e., it does not involve any approximation or truncation with regard to C 0 , in case of Gaussian Eulerian velocities (f c = 0). The solution of Eq. (35) is [17] where g i = g i (u ,x), where a H i = a H i (u ,x) is the solution of the homogeneous problem A variety of solutions exists for a H i , linear and nonlinear in u , but each of them contains a degree of indeterminacy apparent in unspecified constants [17]. When confining the damping function to linear representations in u i , the solution of Eq. (38) is where ε klj is the alternating unit tensor. Solution (39) constitutes an antisymmetric extension to the symmetric damping tensor derived in the previous section and described by the first term of solution (32). In this solution b k are three dimensionless constants whose values are unknown. It is a reflection of the nonuniqueness problem: Except for isotropic turbulence it is impossible to fully specify the damping function on the basis of a specified fixed-point Eulerian velocity distribution [11,15]. Yet there is a practical way out of the nonuniqueness problem. It appears that a H i yields only contributions of relative magnitude O(C 0 −2 ) compared to the previously determined leading terms in the statistical distributions of particle displacement. This conclusion is arrived at when extending the expansion leading to the diffusion equation by one order in C 0 −1 ; see Sec. V below. It reveals contributions of relative magnitude O(C 0 −2 ) in diffusivity and convection only. The same result is obtained for the other term in solution (37) which describes the effect of non-Gaussianity. In general, the contribution of g i in solution (36) can be disregarded in any description which allows for a relative error of O(C 0 −2 ) in dispersion statistics (see the next section). Setting g i = 0 we arrive at a Langevin model which has as a damping function This model can be further reduced to elementary form by replacing the third term by the more simple drift term (∂/∂x k )σ ik , which is the fixed-point statistical average of the third term, From the analysis in the next section, it can be verified that representations (40) and (41) result in the same truncation error of O(C 0 −2 ) in statistical displacement. While the first term corresponds to the result of the Hamiltonian base case, the second and third terms in solution (41) represent the correction due to inhomogeneity in an otherwise locally homogeneous statistical field. The corrections can be related to the change of energy which was disregarded in the leadingorder formulation where underlying particle mechanics were considered Hamiltonian. The second term describes the change of energy due to changes of covariances in the direction of the mean flow. This will happen in turbulence in accelerating or decaying mean flow. The third term describes the effect of the spatial gradient of fluid velocity covariance. This can be associated with shearing due to external forcing.
Solution (40) corresponds to an earlier result of Thomson [4]. It was one of several proposals made for the damping functions which all satisfy the well-mixed criterion and which correspond to an entirely Gaussian Eulerian velocity distribution [13,15]. It is a reflection of indeterminacy because of nonuniqueness. The present analysis provides a solution. It reveals descriptions for statistical displacement obtained from Eqs. (40) or (41) which are unique up to an error of O(C −2 0 ). In Sec. VI we will indicate that this error is small and amounts to only a few percentages for shear-induced turbulence. 066309-6

V. HIGHER-ORDER FORMULATION OF THE DIFFUSION EQUATION
Fokker-Planck equation (30) enables derivation of a diffusion equation for particle dispersion which is one order in C 0 −1 more accurate than Eq. (29). The procedure is similar to that leading to Eq. (29) except that the expansion is extended with one order in C 0 −1 and the full description of a i given by Eqs. (32)- (38) is taken into account. The procedure has been presented in Refs. [16,17]. A different method is given below. Here we start from the conservation equation for passive admixture and derive expressions for mean concentration and mean turbulent flux using the Fokker-Planck equation. The method involves expansion schemes for probability density which are limited to two terms. It is appreciably shorter and more transparent than the solution procedure published previously in which a four-term expansion scheme was employed.
On disregarding transport by molecular motion, a passive scalar χ carried by an incompressible turbulent fluid satisfies the conservation equation The conserved scalar can be fluctuating concentrations of a passive admixture such as smoke or fumes or fluctuating temperatures of an (almost) incompressible fluid [1]. For the Eulerian flow field we can write u i = u 0 i + u i , where u 0 i is mean velocity and u i the fluctuating part of velocity at a fixed point. Now average Eq. (42) at a fixed point; such averaging is indicated by angled brackets. We then have where the mean concentration χ and the mean turbulent flux χu i can be related to the joint probability density p(v ,x) of velocity and position of admixture particle by the equations The density p is described by Fokker-Planck equation (30). In terms of the time t defined by Eq. (17) and a i defined by Eq. (32), the Fokker-Planck equation becomes We now represent p by the expansion of Eq. (18). Substituting this expansion in the above equation and equating terms of leading power in C 0 −1 yields Eq. (19) as description for the leading term p 0 . Its solution is given by Eqs. (21) and (22). For the next-to-leading-order term we obtain, from Eq. (46), Invoking the solution for p 0 , the expressions for a i [cf. Eqs. (36) and Eqs. (37)], and noting that (∂/∂t ) the solution of which is given by where the first term on the right-hand side is the homogeneous solution while the second and third terms are particular solutions. Because p G is zero mean and because of the first of relations (34), integrals +∞ −∞ dv of the particular solutions are zero. Hence, +∞ −∞ p 1 dv = G 1 , so we can write, for the mean concentration defined by Eq. (44), The next step is to derive a two-term description for the righthand side of the equation pertaining to conservation of mean admixture concentration [Eq. (43)]. To that end, an expression is derived for the turbulent flux as defined by Eq. (45) from Fokker-Planck equation (46). Therefore, multiply the left-and right-hand sides of Eq. (46) with v k , integrate with respect to v over its entire domain, apply partial integration to integrals containing derivatives with respect to v i , and use the property that probability densities go exponentially to zero as |v | → ∞. Multiplying the left-and right-hand sides of the resulting equation with σ mk and using σ mk λ kj = δ mj , we obtain, for turbulent flux, the equation where we used the property that integrals of p with respect to v over the entire v domain vary over the slow time t and do not depend on t . Introducing into the right-hand side the expansion for p [cf. Eq. (18)], describing turbulent flux according to and equating terms of equal power in C 0 −1 yields and The integrals on the right-hand sides of the above equations can be evaluated on substituting the expressions for a i , cf. Eqs. (36)- (38), and the solutions for p 0 and p 1 given by Eqs. (21) and (49), and making use of such properties as which are obtained by partial integration, using the property that (∂/∂v i )p G = −λ ij v j p G and implementing relations (34) and (38). It results in where and The terms d mj and e m represent the effect of the "nonunique" solution a H i and of the "non-Gaussianity" function f c in damping function (36) and (37). These terms only contribute in the next-to-leading-order description of turbulent flux.
Equations (50), (52), (56), and (57) constitute two-term descriptions in powers of C 0 −1 for the left-and right-hand sides of the equation of conservation of mean admixture; cf. Eq. (43). Invoking these descriptions into Eq. (43) and equating leading terms on the left-and right-hand sides, respectively, yields Note here that the left-and right-hand sides become of the same order with respect to C 0 −1 once we introduce the slow time t via Eqs. (17) and (23). The above equation is the leading-order description of the diffusion equation. It differs from Eq. (29) by the presence of a single diffusion coefficient 2C 0 −1 ε −1 σ il σ lj . This is due to the presence of the drift term 1 2 λ jn (∂/∂x m )σ ij (u m u n + σ mn ) or, what is equivalent, (∂/∂x k )σ ik , in damping function (36). This term was absent in damping function (15), which was underlying diffusion equation (29). There is no contribution in diffusion equation (60) of the terms related to non-Gaussianity f c and of nonunique solution a H i in damping function (36) and (37). Potential contributions have vanished because of relations (34) and (38). The vanishing of contributions of a H i can also be understood by noting that a H i constitutes an antisymmetric contribution in the correlation function of particle velocity [16]. In case of a homogeneous anisotropic process, i.e., the Hamiltonian base case, diffusion coefficients are determined by the time integral of the combined velocity correlation function according to [1], Antisymmetric contributions therefore cancel out in the leading-order formulation.
Equating terms next to leading order on the left-and right-hand sides of the equation of conservation of the mean admixture yields The tensor γ kn defined by Eq. (58) can be shown to be antisymmetric. This follows from Eq. (38) on multiplying this equation with u m u n and applying partial integration, It implies that d ij is also antisymmetric, d ij + d ji = 0, so The consequence is that the second-to-last term on the righthand side of Eq. (62) reduces to a convection term. The last term is also a convection term. The terms represent the effect of the nonunique solution a H i and of non-Gaussianity f c in the damping function. When we bring these terms to the left-hand side and join them with the other convective terms, they describe relative contributions of O(C 0 −1 ), that is, O(C 0 −2 ) compared to the leading terms of convection. This is beyond the accuracy of the perturbation scheme, which is limited to a two-term expansion for the left-and right-hand sides. Consistency requires dropping these terms. Multiplying the resulting equation with C 0 −1 , adding this to Eq. (60), and combining leading-order and next-to-leadingorder description of the mean concentration by the single concentration G, we arrive at the diffusion equation The above diffusion equation inhibits a relative error due to truncation of O(C 0 −2 ) in contrast to Eq. (60), which is accurate to O(C 0 −1 ). The differences between the two equations are due to the next-to-leading-order contributions in damping function (40). There is no contribution of g i of damping function a i described by Eq. (36) in the above diffusion 066309-8

VI. ACCURACY AND GENERALITY
We derived analytical expressions for damping in the Langevin equation and for diffusivity in the diffusion equation through application of a perturbation expansion in which the inverse of the Kolmogorov constant served as small parameter.
In the subsequent analysis we shall investigate the accuracy of these expressions and their range of applicability. We shall make use of results from theory, experiment, and DNS of actual cases of turbulent flow.

A. Correlations functions of fluid particle velocity
For the perturbation expansion to be effective it is necessary that the first term which corresponds to the Hamiltonian base case is the dominant term in the description of the solution as a whole. This is what is found in the results for autocorrelation functions of particle velocities obtained from theory, measurement, and DNS for three important cases of turbulent flow: decaying grid turbulence, pipe flow, and channel flow. As shown in Figs. 1-3, in all these cases results agree quite well with autocorrelations obtained for the Hamiltonian base case, i.e., the model in which nextto-leading-order terms in the damping function are dropped altogether and where the coefficients of the remaining linear damping term are taken as constant and equal to their values at the point of marking; this is the model according to Eqs. (1), (4), (10), and (14). A derivation of the exact result for decaying isotropic Gaussian turbulence behind a grid can be found in Appendix A. Formulae for the correlation functions appropriate for the Hamiltonian base case are presented in Appendix B. Extensions of these formulas to account for the effect of finite viscosity are given in Appendix C. These extensions are needed in comparisons with results of measurements and DNS at finite Reynolds number (the cases presented in Figs. 2 and 3).
Hamiltonian base case DNS experimental uτRν -1 = 362 2. (Color online) Lagrangian velocity correlations for pipe flow. Dots are experimentally obtained results using particle tracking velocimetry [36], dashes lines are results from DNS [37], and solid lines are results from the Hamiltonian base model corrected for finite viscosity. Input data for the Hamiltonian model are in accordance with DNS data:  [18]. This has been illustrated in Fig. 4, where we have plotted cross-correlation functions of channel flow. Results shown apply to a particular position in the channel; similar agreement has been found at all other distances away from the viscous layer at the wall [18]. The agreement implies that also off-diagonal components of the diffusion tensor can well be predicted using the Hamiltonian base model. Diffusion coefficients determine dispersion. The agreement between velocity correlation functions thus indicates that, for the considered cases, dispersion is already well described by the first term of the expansion, i.e., by a model which corresponds to a locally homogeneous anisotropic Hamiltonian process of velocity fluctuations.

B. Non-Gaussianity and antisymmetry
As mentioned, to improve accuracy we have extended the Hamiltonian base case with next-to-leading-order descriptions. It has resulted in a Langevin model in which Eq. (41) describes the damping function. Application of this model leads to an error of O(C −2 0 ) in statistical displacement. In Eq. (41) we have disregarded the contribution of g i in the general solution of the damping function given by Eqs. (32) and (36)- (38) and obtained through application of the wellmixed principle. Here g i consists of two parts [cf. Eq. (37)], a contribution because of non-Gaussianity described by f c and the contribution a H i which models antisymmetry in the damping function. Values for f c can be inferred from the many data which are available for Eulerian velocity statistics [5,[27][28][29][30][31]. They cover a wide range of cases of turbulent flow and all show mild deviations from Gaussian behavior. Taking values for f c which are in line with these data, one finds corrections in statistical displacement which are only a few percentages [17].
To asses a H i , we can make use of the results of DNS and measurements of turbulent flow in pipes and channels [16][17][18]36,37] and in uniform shear flow [14]. Anisotropy of the Eulerian covariance tensor is a measure for antisymmetry in the damping function. Anisotropy due to shearing is known to be large in pipe and channel flow and uniform shear flow. Values for a H i derived from data obtained for these cases are representative for many, if not all, sorts of turbulence in which shearing plays a role. They can serve as a rigorous test case for the presented results and method of solution when it concerns errors because of disregarding a H i . To calculate the contribution of a H i we can use solution (39), which models antisymmetry. The value of the unknown coefficient b k in this solution can be deduced from the antisymmetry between cross-correlation functions which were measured and calculated by DNS. The effect of the value of b k on dispersion can subsequently be inferred from a diffusion coefficient in which b k is taken into account (the diffusion equation can also be derived when description (39) is included in the damping function; b k then explicitly occurs in the diffusion coefficient [16]). In this way it was found that for the various cases considered the effects of antisymmetry on the diffusion coefficient were minor, i.e., a few percentages at most.
Antisymmetry of limited magnitude was also seen in the results of DNS for channel flow [18]. As illustrated in Fig. 4, the effects of antisymmetry apparent in the difference between the opposing cross-correlation functions were rather small. This behavior was seen at all positions away from the viscous layer at the wall. Mild forms of non-Gaussianity and antisymmetry comply with the outcome of the present analysis. Energy change takes place on a time scale which is appreciably longer than that of statistical velocity fluctuation (by a factor C 0 ) so underlying particle mechanics are approximately Hamiltonian

C. Analytical evidence
Statistical turbulence is known to be a basically inhomogeneous process in which there is continuous supply and dissipation of energy. When following a marked fluid particle, statistical averages of fluid velocity and energy change significantly in magnitude over the time scale σ ε −1 [e.g., cf. Eq. (8)], where σ is magnitude of mean-square fluctuating velocity. As follows from Eqs. (1) and (15) to distinguish a shorter time scale where statistical velocity can be modeled according to a homogeneous anisotropic Hamiltonian-based process. That this is possible can most directly be inferred from the expressions in closed-form for decaying grid turbulence presented in Appendix A. Here it can be seen that expansion of the exact results in terms of powers of C 0 −1 leads to the Hamiltonian base case as the leading-order description. As shown in Fig. 1, the leading-order description is already close to the exact result.
For pipe and channel flow it also can be shown in more explicit analytical terms that approximation by expansion in powers of C 0 −1 works [16]. For wall-induced turbulence there exists an asymptotically exact description for ε in the log layer [38]. The correctness of this expression has, among others, been confirmed by measurements in the atmospheric surface layer [39]. Although in the log layer ε is inversely proportional to distance from the wall, it changes only C 0 −1 on the time scale of the particle velocity correlation [16]. It is for this reason, among others, that, even in the log layer where ε changes significantly in magnitude, the correlation functions obtained from DNS are well approximated by those of the Hamiltonian base model in which ε is kept constant and equal to its value at the point of marking; i.e., the model according to Eqs. (1), (4), (10), and (11). Moreover, the approximation of disregarding antisymmetry in the damping function when describing dispersion seems to be correct. All this is in line with an approximation scheme based on small C 0 −1 and where the Hamiltonian base model arrives as a leading-order formulation from the general model of Eqs. (1), (32), and (36)- (38).

D. Higher-order terms and dispersion
The next-to-leading-order terms in the damping function which are to be retained in a two-term expansion for statistical displacement are the second and third terms on the right-hand side of Eq. (41). Both terms describe the effect of covariance inhomogeneity on the short time scale of particle velocity fluctuation. The last term results in mean particle drift, and the second term leads to a correction in diffusivity as can be seen from the second term in solution (67). Although the mean flow u 0 n is generally large in turbulent flow, the magnitude of the second term is limited because of conservation of momentum. An illustration of this is decaying turbulence behind a grid, for which exact results are available [17]. For grid turbulence the second term in Eq. (67) is 2 3 C −1 0 times the first term. Comparing the two-term description for diffusivity according to Eq. (67) with the exact result, we find a relative deviation of 4 9 C −2 0 which amounts to 1.2% when C 0 = 6. It shows that the error due to truncation of O(C −2 0 ) is limited in the diffusion equation as well. An exception is the log layer of wall-induced turbulence where the error entailed by the diffusion model is somewhat larger (about 10%). This is because of the strong variation of ε and the diffusion coefficient with distance from the wall [16].

E. General shear-induced turbulent flow
The two-term descriptions in powers of C 0 −1 constitute a description for statistical particle displacement which inhibit a relative error due to truncation of O(C 0 −2 ). For the above considered cases of decaying grid turbulence, uniform shear flow, pipe flow, and channel flow this error was found to amount to a few percentages. The error does not differ much from the value of C 0 −2 itself. In view of the general principles which underly the derivation, the truncation error will also be O(C 0 −2 ) in other cases of turbulence. More specifically, we can expect errors of a few percentages in all those cases where the conditions in terms of anisotropy, inhomogeneity, non-Gaussianity, and the value of C 0 −1 are similar to those of the considered cases. This will be true for many, if not all, sorts of turbulence which originate from shearing. In addition to the cases considered we mention turbulent flow in boundary layers, including the atmospheric surface layer, turbulent jets and wakes, Couette flow, and flow between counter-rotating disks. In all these cases, reduction of the general Langevin model of Eqs. (1), (32), and (36)- (38) to the model of Eqs. (1) and (41) is expected to lead to limited error in the prediction of displacement statistics. The linear part of the damping function is the important term under all circumstances. This is because of its functional form. It contains the coefficients C 0 , ε, and λ ij , where C 0 is relatively large while ε and λ ij always have appreciable values whatever the sort of turbulence. All other terms are less significant in magnitude. The contributions due to non-Gaussianity, antisymmetry, and inhomogeneity are expected to not differ much from those identified for wall-induced turbulence in pipes and channels where they were found to be small. Moreover, the inverse of the Kolmogorov constant can be used as a parameter to measure and order terms in the statistical model.

VII. CONCLUDING REMARKS
In our derivation of the statistical description of turbulent dispersion we started from the Markov approximation for particle velocity. It finds its justification in Lagrangian Kolmogorov theory, the outcome of which is supported by results of DNS and measurements [7][8][9]. Given a Markovian-based Langevin equation for statistical particle velocity, at issue is the specification of the damping function. The second issue is the description of statistical particle displacement through a diffusion equation. Both issues were dealt with by employing a perturbation scheme based on expansion in powers of C 0 −1 where C 0 is the universal Kolmogorov constant. Because of limited smallness of C 0 −1 , viz. C 0 −1 ≈ 1 6 , we did not restrict our analysis to terms of leading order but also considered terms next to leading order. The final results are the Langevin model of Eqs. (1), (2), and (41) and the diffusion equation according to Eqs. (66) and (67). They inhibit a truncation error of relative magnitude O(C 0 −2 ) in statistical displacement. This was shown to lead to a deviation of a few percentages in predicted dispersion statistics for representative cases of turbulent flow. The descriptions do not contain fitting constants or empirically established functions. All coefficients are uniquely determined by fixed-point Eulerian-based statistical averages of the flow field and the universal Kolmogorov constant. The models satisfy the condition of being falsifiable.
Perhaps surprisingly, the leading order formulation of the Langevin model was found to correspond to a locally homogeneous anisotropic Hamiltonian-based statistical model of fluid particle velocity. It substantiates the proposal of the 066309-11 pioneers in turbulence theory (Von Kármán, Taylor, Prandtl) to model turbulent dispersion analogous to diffusion by molecular motion [1]. What is new is that the diffusion coefficient is now a well-defined tensor, cf. Eq. (67), and that the accuracy of the model is known. It is noted that the present analysis does not entail a conclusion regarding validity of the concept of turbulent or eddy viscosity.
Allowing for an error of a few percentages in statistical displacement, the presented models can be widely applied to cases of shear-induced turbulent flow. Potential areas for application are diverse and can be of great societal importance. Examples include the following: fumes from chimneys, radioactive aerosols from nuclear accidents, sand particles from dust storms, micron-sized particles from volcanic eruptions, oil spills in the sea, smells from factories, spray injection in technical devices, and so on. The presented descriptions can be applied once the fixed-point Eulerian mean and covariance and the energy dissipation rate of the flow field are known. Such information is available or can be made available through measurements, known analytical descriptions, and results of computer calculations. Dispersion statistics can then be calculated from the Langevin model given by Eqs. (1), (2), and (41) and diffusion equation (66) with diffusion coefficient (67). The Langevin model yields the most accurate results but requires the largest calculation effort. The diffusion equation is easier to handle but is somewhat less accurate.