Statistics of a passive scalar in homogeneous turbulence

Statistics of a passive scalar with Sc=1 transported by steady homogeneous turbulence at Rλ=427 and Pλ=427 is studied by using high-resolution direct numerical simulation. The Obukhov–Corrsin constant of the three-dimensional scalar spectrum in the inertial-convective range is found to be 0.68±0.04. It is proved that the -law for the scalar-velocity triple correlation holds in both inertial-convective and viscous-convective ranges when Sc>1, and found that the -law is approached with increase in Péclet number. Structure functions of the passive scalar increment and their local scaling exponents are computed as functions of the separation distance, and it is found that there exist two scaling ranges: the inertial-convective range and a narrow precursory range to the viscous-convective range. The scaling exponents in the inertial-convective range are found to be smaller than those of the velocity field and do not saturate, whereas they saturate at about 1.5 in the short precursory range to the viscous-convective range. It is also found that, contrary to the scalar case, the mixed scalar velocity structure function has a well-developed single scaling range. The scalar and scalar dissipation fields are visualized and compared with the kinetic energy dissipation field. The scalar field has a particular shape with a large-scale plateau, sharp cliff and deep valley, a mesa-canyon structure.


Introduction
The Kolmogorov theory is a milestone for the statistical theory of turbulence [1,2]. With two hypotheses, Kolmogorov obtained δu r ∝ (¯ r) 2/3 scaling for the velocity increment δu r = u(x + re 1 ) − u(x) in the inertial range η r L u , where L u is a macroscale, η = (ν 3 /¯ ) 1/4 is the Kolmogorov scale, ν the kinematic viscosity, and¯ is the mean rate of the energy dissipation of a unit mass of fluid. Then the energy spectrum is given by where K is a non-dimensional constant and has been estimated to be ∼1.62 by experiments and 1.64 by direct numerical simulation (DNS) [3,4]. Kolmogorov also derived the 4 5 -law, which is asymptotically exact and is now used to define the inertial range in data obtained from experiments and DNS. Passive scalar transport by incompressible turbulence has been attracting considerable interest owing to its important applications in industrial machine design and also due to its energy and contaminant transport in geophysical and environmental research [2,5]. By applying Kolmogorov's ideas to the passive scalar, Obukhov [6] and Corrsin [7] derived the scaling law δθ 2 r ∝χ¯ −1/3 r 2/3 (3) for the second-order moments of the scalar increment δθ r = θ(x + r) − θ(x). The scalar spectrum is defined by and the scalar spectrum in the inertial-convective range η r L θ is given by where L θ is the integral scale of the scalar andχ is the average rate of smearing of the passive scalar by molecular diffusion; we call it the scalar dissipation rate throughout this paper. The non-dimensional constant C OC is the Obukhov-Corrsin constant and considered to be universal. Sreenivasan [8,9] examined carefully all the passive scalar spectra in the literature and argued that, in the grid turbulence, the inertial-convective range spectrum of the scalar can be discussed even at modest Reynolds numbers, while in the shear flow turbulence the scaling exponent of the scalar spectrum E θ ∝ k −m θ depends on the Reynolds number and approaches 5 3 when R λ > 1000. This suggests that the homogeneous isotropic scalar turbulence reaches an asymptotic state suggested by Obukhov and Corrsin at a rate faster than in the shear-flow turbulence; Sreenivasan concluded that the Obukhov-Corrsin constant is about C OC = 0.68 [2], [8]- [10].
When the Schmidt number Sc = ν/κ > 1, the energy spectrum decays quickly at wavenumbers larger than k d = 1/η, whereas the scalar spectrum remains excited at levels higher than the energy. The scalar transfer to small scales is maintained by the action of velocity straining whose characteristic time is (ν/¯ ) 1/2 , independent of scale, and is ceased at the scale η B = (κ 2 ν/¯ ) 1/4 by molecular diffusion. This range of scales is the viscous-convective range and the scalar spectrum obeys a k −1 power law where B is a non-dimensional constant [11,12]. Non-locality of the scalar transfer in the wavenumber space is essential for the generation of the viscous-convective range. The value B is also considered to be universal and was about 3-6 [13]- [21]. A relation analogous to the 4 5 -law for the passive scalar was derived by Yaglom [22], the 4 3 -law, δu r δθ 2 r = − 4 3χ r for η B r L θ .
It is important to notice that the lower end of the range over which the 4 3 -law holds extends up to η B ( η) when Sc 1 because η B = Sc −1/2 η. This means that although the passive scalar spectrum has a power (− 5 3 ) in the inertial-convective range and (−1) in the viscous-convective range, the scalar flux in the wavenumber space θ (k), which is defined as the flux of the scalar variance across k transferred to wavenumbers higher than k, is independent of the wavenumber in both ranges. In other words, the scalar flux has a single scaling range η B r L θ which is wider than the case of the scalar alone.
One of the most important issues in turbulence research is intermittency. When the Reynolds number is sufficiently large, the intermittency of the velocity field builds up with decrease of scale, and the probability distribution of the velocity field changes shape. If Kolmogorov scaling holds, the scaling exponent of the structure function for the velocity increments defined by (δu) q ∝ r ζ q is given by ζ q = q/3, but experimental and DNS data have shown that ζ q is a non-decreasing function of q rather than a linear function of q [23]- [31].
The passive scalar transported by turbulence also exhibits intermittency, and its strong fluctuations are in general related to important and/or serious practical problems such as mixing in chemical reactions and localized high concentration of air pollutants and so on. Studies have shown that the scaling exponent of the passive scalar is smaller than that of the velocity, meaning that the intermittency of the scalar field is stronger than the velocity field. There has been considerable progress in understanding the statistical properties of the Kraichnan model for the passive scalar, where the scalar is convected by a Gaussian random-velocity field with very short correlation time [13,14], [32]- [37]. Anomalous scaling exponents have been derived analytically for some limiting cases and saturation of the scaling exponents at higher order is suggested [38,39]. However, since the velocity field in the Kraichnan model is artificial, it is uncertain whether a similar behaviour of the anomalous scaling exponents holds for the case of the passive scalar convected by a generic Navier-Stokes turbulent flow [5], [40]- [42]. Furthermore, the importance of ramp-cliff structures and anisotropy effects on intermittency have been recognized when a mean scalar gradient exists. Ramp-cliff structures effectively produce large scalar gradients with strong scalar dissipation, and such structures induce anisotropy lingering at scales below which the velocity field remains isotropic. These factors enhance intermittency. It is, however, not known whether such ramp-cliff structures exist when a mean uniform scalar gradient is absent.
High-resolution DNS of the passive scalar provides us with a very useful method to study various aspects of the scalar problem. Regarding DNS study for the passive scalar problem, it is important to note the following facts. Since Pe = Sc Re, an asymptotic state Pe 1 implies three regimes: regime 1, Re = O(1) or moderate and Sc 1; regime 2, Re 1 and Sc = O(1) or moderate; regime 3, Re 1 and Pe 1. Regime 1 has been extensively analysed in detail, especially with emphasis on the scalar spectrum and intermittency of the scalar gradients, because there is almost no 'bump' in the passive scalar spectrum in the crossover region between the viscous-convective and diffusive ranges, and most of the wavenumber range in DNS for the passive scalar is effectively assigned to resolve the viscous-convective range; thus a welldeveloped scalar spectrum E θ ∝ k −1 has been obtained easily when compared with regime 2. In regime 1, several authors have computed the Batchelor constant B of equation (6) and examined the asymptotic form of E θ (k) in the far diffusive range [13]- [21]. The Batchelor constant B is found to be between 3 and 6, and Kraichnan's spectral form E θ (k) ∝ exp(−ckη B ) is found in three dimensions, where c is a non-dimensional constant. Schumacher et al [43,44] studied analytically the Schmidt-number dependence of derivative moments under a mean scalar gradient and found that the upper bound for the nth-order moments of the scalar gradient grows as Sc n/2 at large Sc. They also found an increasing tendency to isotropy with increase in Sc [20].
Regime 2 has been studied by many experimentalists [2,5,8] and recent high-resolution DNS results were quite interesting. The central feature here is the intermittency of the scalar increments in the inertial-convective range. The resolution requirement is restrictive in this case, and the number of DNS studies in this regime is limited. Chen and Cao [45] have computed the scaling exponents of the passive scalar by using DNS with resolution up to 512 3 and found that they are smaller than the velocity field. However, since the width of the inertial-convective range is not wide enough to observe a clear scaling behaviour, the scaling exponents computed must be understood carefully. Wang et al [46] studied the passive scalar statistics from the viewpoint of the refined self-similarity with resolution N = 512 3 up to R λ = 195 and Sc = 0.7 and 1, but the actual time span for the time average was very short. They found that C OC = 0.87 and the intermittency exponent 0.43-0.77. Log-normal distribution of the spatially averaged scalardissipation rate was observed. Yeung et al [20] found that C OC is ∼0.67 at R λ = 240 and Sc = 1 under a mean scalar gradient. However, as stated below equation (6), when the non-locality of the scalar transfer in the scale space is strong, there are possibilities that anisotropy at large scales remains at small scales and that the lower end of the inertial-convective range is masked by the excitation of the scalar amplitudes in the viscous-convective range prior to roll-off of the scalar spectrum due to the molecular diffusivity. This occurs when the resolution range in DNS is not sufficient, and may lead to a spurious wider inertial-convective range, so that the scaling exponent in the inertial-convective range thus computed may suffer from systematic deviation. To avoid such possibilities, it is desirable to compute the structure functions with sufficient width of scale range such that the lower end of the inertial convective range is clearly discriminated from other scaling ranges when Sc 1. In this context, the present DNS study is the first, and is capable of analysing the scaling behaviour of the passive scalar in detail. Also this DNS provides data which enable us to determine the scalar flux transferred from the inertialconvective range to scales smaller than that.
Regime 3 is far beyond the scope of DNS by present day computers including Earth Simulator [47]. Only an experimental approach is available.
Our purpose here is to present new statistical data on the passive scalar by using veryhigh-resolution DNS, especially focusing on the small scales, and to analyse them from the viewpoint of their scaling behaviour and intermittency with special emphasis on the inertialconvective range. We explore regime 2 and the reason is two-fold. First, in the comparison of the velocity and scalar statistics, Sc = 1 means that the contributions of the molecular dissipation and diffusion to the scalar transfer in the scale space and to the intermittency become equal. Then the difference between the two statistics can be interpreted as due to the pressure term and to the incompressibility of the velocity field. We can examine the effects of the Schmidt number on the various statistics of the scalar by comparing with the case of Sc = 1 as a reference.
Secondly, in the literature, there are no systematic DNS data for the structure functions, local scaling exponents and probability distribution functions (PDFs) of scalar increments in the inertial-convective range. We have already made a high-Reynolds-number DNS, R λ = 460 [4]. The results and experience of the velocity statistics by that DNS are used effectively in many ways to study regime 2. The present DNS provides such data in detail. We analyse them carefully, compare their statistics with those of the velocity and try to understand the properties of the intermittency of the passive scalar convected by turbulence. Also, the statistics of the scalar variance flux are examined, and the fields of the scalar and scalar dissipation are visualized and their structure examined.

Numerical simulation
The Navier-Stokes and continuity equations, and the passive scalar equation were numerically integrated under the periodic boundary condition with periodicity 2π. The pseudo-spectral method was used for the non-linear and convective terms, and the time integration was done by the fourth-order Runge-Kutta-Gill method. The Gaussian random solenoidal force f and scalar source g θ are delta-correlated in time and added in the low wavenumber range, 1 k 2. The Schmidt number is unity. Spectra of kinetic energy and scalar variance are defined by The integral scales are defined by the average energy dissipation rate and scalar dissipation rate bȳ the skewness for the velocity and scalar by the Taylor microscale for the velocity and the scalar by and microscale Reynolds and Péclet numbers by The spatial resolutions are N = 512 3 for run 1 and N = 1024 3 for run 2. R λ = P λ = 258 for run 1 and 427 for run 2, respectively. The initial conditions for the velocity and scalar are Gaussian random fields which have the same initial spectra E(k, 0) = E θ (k, 0) ∝ k 4 exp(−2(k/k p ) 2 ) at low resolution N = 256 3 but are statistically independent of each other. After reaching a steady state, (R λ , P λ ) and the resolution were increased to R λ = P λ = 258 with N = 512 3 for run 1. For run 2, we used the velocity field with R λ = 460 of Gotoh et al [4] but with a slightly bigger viscosity and the scalar field of run 1 as the initial condition. This reduced the computational cost. After a transient, a statistically steady state was obtained. Statistical averages were taken as time averages during the steady state, T av = 6T eddy for run 1 and 2.5T eddy for run 2, where T eddy = L u /u rms . The numerical parameters are listed in table 1. A measure of the accuracy of the present DNS is given by K max η = K max η B ≈ 1 which is slightly smaller than that commonly used in other studies [18,20,46]. The spatial resolution of DNS will be discussed in section 6.

One-point statistics
Time evolution of the total energy, total variance, mean rates of the kinetic energy and scalar dissipation were monitored throughout the computation and used to define the steady state. Figure 1 shows the time evolution of the normalized energy and scalar dissipation,ˆ =¯ L u /u 3 rms andχ =χL u /(u rms θ 2 rms ). It is seen that a statistically steady state is established after t = 6 for run 1 and t = 1 for run 2.
As in the case of the energy dissipation, it is expected that, when P λ → ∞,χ remains finite [47]. The normalized valueχ =χL u /(u rms θ 2 rms ) for various R λ is compared with the data from DNS [46] and experiments [10,48] in figure 2. Note that R λ = P λ because Sc = 1. It is interesting to see that, although the number of data points is small, there appear to be two groups of values ofχ for Sc = 1 and 0.7, and each stays approximately at the same constant level. The average value is 0.318 for Sc = 1 and 0.472 for Sc = 0.7, respectively. The ratio of the former to the latter, 0.318/0.472 = 0.67 is very close to 0.7. It is important to examine whether χ tends to an asymptotic value. (t),   Figure 3 shows the compensated three-dimensional spectra of the kinetic energy and scalar variance normalized using Kolmogorov-Obukhov-Corrsin variables. Note that η = η B because Sc = 1. Both E(k) and E θ (k) increase in the range kη 1, implying that those scales are not well resolved. This is due to the finiteness of the wavenumber domain, but the statistics of the two fields at wavenumbers below kη = 1 do not suffer significant contamination of numerical

Spectra
The value K = 1.61 is consistent with the experimental value K exp = 1.62 and the DNS value K DNS = 1.64 [3,4], and C OC = 0.68 agrees well with the recommended value 0.68 by Sreenivasan who examined carefully the values reported from many experiments [9]. The value C OC = 0.68 is also consistent with C OC = 0.75-0.92 (Sc = 0.7, R λ = 582, P λ = 407) obtained from wind-tunnel experiments with uniform temperature gradient by Mydlarski and Warhaft [10]. When compared with other DNS results, Wang et al [46] performed DNS with lowband forcing and found C OC = 0.75 (Sc = 1, R λ = P λ = 195), and Yeung et al [20] obtained 0.67 (Sc = 1, R λ = P λ = 240) with uniform mean scalar gradient. However, in the calculation of Wang et al, the range of forcing overlaps the − 5 3 range, and in the calculation of Yeung et al, the uniform mean gradient may induce some degree of anisotropy of the scalar field in the inertial-convective range. However, in the present DNS, the k −5/3 spectrum exists for wavenumbers higher than the forcing range. We conclude that the horizontal portion of the curves in the present DNS certainly shows the inertial-convective range spectrum of the passive scalar. When compared with the curves for E(k), the curves for E θ (k) have large bumps at kη 0.2, and their peaks shift towards higher wavenumbers compared with the velocity case. This implies that, even when Sc = 1, the non-locality of the scalar variance transfer is stronger than in the velocity case, and the large bump is understood as a precursor of a k −1 spectrum in the viscous-convective range for the passive scalar [49].   [46]. Figure 4 shows comparison of the energy and scalar dissipation spectra 2νk 2 E(k)/¯ η and 2κk 2 E θ (k)/χη B normalized by the Kolmogorov variables. The curves are smooth and collapse well at wavenumbers below kη < 0.8, and the peak points of the spectra are about 2.3 at kη ≈ 0.16 for the velocity and about 1.56 at kη ≈ 0.25 for the scalar, which are the same as those found in figure 2 of Wang et al [46] (note that our definition of D(k) is twice as large as that of Wang et al). On the other hand, the curves at kη > 0.8 slightly rise and do not collapse, suggesting that those scales are not resolved properly.
The energy and scalar variance transfer functions are defined by where F(k) and F θ (k) are the spectra of the random force and scalar source, respectively. Transfer fluxes of the energy and the scalar variance across the wavenumber k are defined by  over which E(k) and E θ (k) both have a slope of − 5 3 , and the constancy of the curves for those wavenumbers means that the turbulent and scalar fields are really in Kolmogorov's equilibrium state. Note that flat regions of θ curves extend to wavenumbers higher than in the velocity case.

Second-order structure functions
The structure functions for longitudinal velocity increments δu r , transverse velocity increments δv r and the scalar increment δθ r are defined by In the Kolmogorov theory, these three structure functions scale as r q/3 . Figure 6 shows the second-order structure functions¯ −2/3 r −2/3 S L,T 2 (r) andχ −1¯ 1/3 r −2/3 S θ 2 (r). As found in [4], the scaling exponent of the velocity at the second order is about 0.696, which is larger than 2 3 , and the curves for the velocity compensated by multiplying by r 2/3 are slightly increasing as r increases. The wavenumber range 0.008 kη 0.03 over which E θ (k) ∝ k −5/3 corresponds to the scale range 209 r/η 785, and in fact we see that the compensated curve for S θ 2 is almost horizontal. When r decreases, the curve for the scalar in run 2 changes slope slowly over the range 100 r/η 200 and reaches a peak and then quickly decreases. The crossover range 100 r/η 200 corresponds to the wavenumber range 0.03 < kη < 0.06 over which the crossover of E θ (k) to the bump occurs. These observations show that there is a oneto-one correspondence between the spectrum and the second-order structure function of the passive scalar.
The isotropy of the velocity field in the present type of DNS was well established at the level of the second-and third-order moments at scales below about half the integral scale L u [4]. 10 100 1000 For example, examination of the longitudinal and transverse structure functions reveals that they satisfy the relation S T 2 (r) = S L 2 (r) + (r/2)∂S L 2 (r)/∂r (figure not shown), which is derived under the constraints of incompressibility and isotropy. Since the incompressibility condition is well satisfied in our DNS to the order of machine precision, we conclude that the isotropy of the velocity field is also well achieved, as in the simulations of Gotoh et al [4].
The local scaling exponents of the structure functions are defined by where α denotes L, T or θ.  r/η = 300. Although it is not certain whether the flat portion in the ζ θ 2 (r) curve emerges at the level of those local minimum and maximum values as Péclet number increases, the local maximum value 0.66 is already within the 5% error bound of the velocity scaling exponent 0.696. If the Kolmogorov theory does not apply to the second-order scaling exponents, there seems no reason to expect ζ θ 2 = ζ L 2 = ζ T 2 . The local maximum and minimum in ζ θ 2 (r) imply the existence of two scaling ranges. 4 3 -law for the scalar variance Kolmogorov's 4 5 -law is an asymptotically exact result for turbulence, and the corresponding formula for the passive scalar was derived by Yaglom [22]. When the Reynolds number and Péclet number are very large, the scalar field in a steady state obeys the equation

3.4.
where F θ (r) is the input term from the random scalar source acting at large scale. When the separation distance r is much smaller than the macroscale L θ and much larger than the diffusive scale η B , the second and third terms in (22) are negligible, thus the 4 3 -law is obtained as It is important to notice that the lower end of the range for the 4 3 -law to hold is not at η but at η B . Since η B = Sc −1/2 η for Sc 1, the above 4 3 -law is valid for both inertial-convective and viscous-convective ranges when Sc 1. Figure 8 presents curves of − δu 3 r /(¯ r) for the 4 5 -law and − δθ 2 r δu r /(χr) for the 4 3 -law [4,10,20,48,50,51]. When R λ and P λ increase, the width of the plateau of the curves becomes wider and their values approach 4 5 and 4 3 , respectively. The plateau of the curves for the scalar extends to scales slightly smaller than those of the velocity field, which is consistent with theoretical expectations and with observation in figure 5. This can be compared with the results of Yeung et al [20] since the 4 3 -law is approached when Sc increases from 1 to 64 with R λ = 38. These observations are consistent with the above statement about the effective range of the 4 3 -law. To see the degree of the fluctuations in the curves, it is useful to compare the curves for − δu r δθ 2 r /(χr) averaged over five successive time spans T i = 0.5T eddy (i = 1, . . . , 5) in figure 9. If we take the curve for total average T = 5T eddy as a reference, then we find that the fluctuation level is about 10% and stronger than in the velocity case. Since the third-order moment is given by the difference between the positive and negative contributions of the PDF for δu r δθ 2 r , where the sign change is due to change of sign of δu r , the fluctuation level is larger than it is in a positive-definite case such as δθ 2 r .

Probability distribution function
The intermittency in hydrodynamic and scalar turbulence is characterized by the growth of the tail of the PDF of the velocity and scalar increments. Figure 10 shows the one-point PDFs for the velocity in the x-direction and the passive scalar. The PDF P(u x ) is very nearly Gaussian. For normalized scalar amplitudes smaller than 4, P(θ) is also very close to Gaussian, but for larger amplitudes it decays smoothly and faster than the velocity PDF. From our experience with DNS for the scalar we feel that in order to obtain a reliable one-point PDF for the scalar amplitude, averages must be taken over a sufficiently long time period, and we must be careful to resolve the tail of P(θ). The asymptotic form of the PDF for large scalar amplitudes is difficult to determine. The one-point PDFs for the gradient fields of the velocity and scalar are shown in figure 11. The PDF for the longitudinal velocity gradient ∂u/∂x is negatively skewed, whereas that for the transverse velocity gradient ∂u/∂y is symmetric and has a tail longer than in the longitudinal case. The PDF for the scalar gradient is higher than those of P(∂u/∂x)   and P(∂u/∂y) for the normalized amplitudes between 5 and 25, but the far tail region decays faster than the velocity gradient PDFs.
To see how the distributions of the velocity and scalar field change with decrease in scale, we plot in figures 12-14 the PDFs for the longitudinal velocity increment δu r , transverse velocity increment δv r and the scalar increment δθ r for run 2 at the separation distances r n /η = 2 n−1 (dx/η), n = 1, 2, . . . , 10, where dx = 2π/N.
When the separation distance r becomes smaller, the tails of the PDFs become longer; P(δθ r , r) is almost symmetric. The PDF tails for the velocity increments are concave, whereas those for the scalar increments become convex at large amplitudes. This feature of the PDF of the scalar is consistent with the experimental and DNS observations, but our PDF    curves continue smoothly to the order of 10 −10 . The scaling exponents of δθ r in the range 30 < r/η < 60 are seen to be approximately saturated as q becomes large as discussed in section 5.2 [41,42]. This suggests that the tail of P(δθ r , r) in this range is of the form of P(δθ r , r) = θ −1 rms r η ζ θ,VCP ∞ Q δθ r θ rms for δθ r δθ r 2 1/2 ,    where Q(x) is a non-dimensional function and ζ θ,VCP ∞ ≈ 1.5 is a value of ζ θ q extrapolated to large q. The scaled PDF curves are plotted in figure 15, where r n /η = 2 n−1 dx/η = 22.3, 44.6, 89.2 for n = 4, 5, 6 are well within the range 20 < r/η < 90. Collapse of the PDF tails is very satisfactory, which means that the large scalar increment scales as θ rms in this range.

Skewness and flatness
The change of the distributions of the velocity and scalar fields with decrease in scale is also characterized by the structure functions. The skewness and flatness of the scalar increment defined by and those of the velocity increments are plotted in figures 16 and 17. When the separation r/L α decreases, the skewness of δu L r becomes more negative, which is consistent with the  5 -law, but that of the transverse velocity is zero. The skewness of the scalar for run 1 is nearly zero, while it is a small positive constant at all scales for run 2. This is due to the large-scale anisotropy, and when a longer time average is taken, the skewness vanishes. This can be compared with the fact that when there is a uniform scalar gradient, the skewness of the scalar increment parallel to the mean scalar gradient is positive and increases with decrease of r [10,20,43,44]. As for the flatness, all the curves increase from 3, the Gaussian value, when r/L α decreases, meaning that the intermittency of the scalar increment becomes stronger. The value of the scalar flatness and its rate of increase are greater than those of the velocities. The value of the scalar flatness for run 2 at the smallest value of r/L θ ≈ 5 × 10 −3 (which corresponds to r/η ≈ 3) is 20, and close to the experimental value of ≈20 [10] and to the DNS value 21.2 for the scalar difference along the mean scalar gradient [20]. Figures 18 and 19 show S L q (r), and ζ L q (r), and figures 20 and 21 show S T q (r) and ζ T q (r), respectively. As found by Gotoh et al [4], there are flat regions of finite width in the scaling exponents for the velocity. On the other hand, figures 22 and 23 show S θ q (r) and ζ θ q (r), respectively, from which there appear to be two scaling regions in which the local scalar scaling exponent takes a local minimum and a local maximum. The crossover occurs at about 100 < r/η < 200. The difference between the local maximum and minimum in the range 20 r/η 600 is small when the order of the moments is small, but gradually increases with the order. In the scaling range 300 r/η 600, the local maximum values of ζ θ q (r) increases with the order but the rate of increase becomes smaller, while for 20 < r/η < 90, ζ θ q (r) (or local minimum values) approaches an asymptotic value about ζ θ,VCP ∞ ≈ 1.5. The former range is the inertial-convective range. Although Sc = 1 is too low for the viscous-convective range to appear, our observation in figures 22 and 23 indicates that the scaling behaviour of S θ q (r) in the range 20 < r/η < 90 is  different from the one in the inertial-convective range, and that S θ q (r) in this range has values higher than those extrapolated from the lower end of the inertial-convective range to smaller scales. These facts are interpreted as the non-local effects of scalar transfer in the scale space, which are the essence of the dynamics in the viscous-convective range, and suggest that the nonlocal effects of the passive scaler transfer are more significant than we would expect and they are appreciable even when Sc = 1. For this reason, we call this short range the viscous-convective precursor range in this paper. It is highly probable that macroscale conditions such as anisotropy  will affect the small scales of the scalar even when the macroscale effects on the velocity field already vanish.

Structure functions at high order
Although it is difficult to draw definite conclusions from these computations, it would be safe to state that the scaling exponents of the passive scalar so far obtained in both DNS and experiments should be examined with great care. Especially when Sc 1 and the width of the scaling range of the structure functions is not long enough, one might mistake ζ θ q in the viscous-convective range for its value in the inertial-convective range, so that values at high q would either be underestimated or misread as showing that the scaling exponent in the inertial convective range becomes saturated.  of 60 < r/η < 400 at low order and its width gradually deceases with increase in q. In other words, it appears that for a fixed q, when Sc 1, S θL q (r) scales with a single power regardless of whether r lies in the inertial-convective or the viscous-convective precursor range. It is, however, not certain whether such scaling behaviour of S θL q (r) continues to prevail when P λ and Sc become very large.
As discussed in the subsection of the 4 3 -law, δu r δθ 2 r = − 4 3χ r holds in both the inertialconvective and viscous-convective ranges. This reflects the fact that the scalar flux θ (r) ∼ δu r δθ 2 r /r is independent of r and equal toχ in both ranges. If we estimate δu r δθ 2 r in the Kolmogorov theory, we have δu r ∼ (¯ r) 1/3 and (δθ r ) 2 ∼¯ −1/3χ r 2/3 in the inertial-convective range, giving δu r δθ 2 r ∝χr, while in the viscous-convective range we have δu r ∼ (¯ /ν) 1/2 r, so that δθ 2 r ∼ (¯ /ν) −1/2χ is independent of r. This is consistent with the k −1 spectrum of the passive scalar because the Fourier transform yields δθ 2 r ∝ ln(r/η) and the logarithmic dependence can be treated as constant for r/η 1. The decay of the velocity amplitude in the far dissipation range is compensated by the constancy of the scalar increment δθ 2 (r = η) for r in the viscous-convective range. The above observation motivates the following conjecture: the scalar transfer function at a given order has a single scaling exponent throughout the inertialconvective and viscous-convective ranges. It seems less plausible that this conjecture would be true at all orders, but it is possible that it holds approximately for low orders. Critical examination by theory, DNS and experiments is certainly required to evaluate this conjecture. Figure 26 presents ζ α 6 for α = L, T, θ, θL and θT . The curves except ζ θ 6 are almost constant in the range 100 < r/η < 300 which corresponds to the inertial-convective range, and their values are found to be ζ L 6 = 1.75 ± 0.01, ζ θL 6 = 1.55 ± 0.01, ζ θT 6 = 1.54 ± 0.01, (27) from which we can compute the intermittency exponents as The value for the velocity 0.25 is the same as the experimental value. On the other hand, the value for the scalar 0.45 is within the range µ θ = 0.43-0.77 by DNS of Wang et al [46,66] where µ θ was defined by χ 2 r ∝ r µ θ and χ r is the volume-averaged scalar dissipation rate, and also within 0.35-0.72 surveyed by Sreenivasan et al [52]. The recent estimate for µ θ is about 0.35 [53,54,64].

Reynolds number and Péclet number effects
To see how fast the local scaling exponents approach asymptotic values when R λ and P λ are increased, we plotted in figures 27 and 28 the variation of the local scaling exponents for runs 1 and 2 (for those of the velocity, readers may see figures 21-24 of Gotoh et al [4]). Both the local minimum and maximum values of the ζ θ q curves for a fixed q tend to be smaller and their positions shift towards larger separation distance. The curves for ζ θL q extend their horizontal portions towards larger r as (R λ , P λ ) increases, but the flat portions stay at the same values. The rate at which the local maxima become flat seems to be very slow when compared with the velocity case.

Scaling exponents
The scaling exponents are determined by taking averages of the values of the local scaling exponent curves which are at a plateau. However, it is difficult to determine the scaling exponents for the passive scalar, and we consider that it is not appropriate to evaluate them by averaging over the single range of scales 20 < r/η < 600. Instead, we computed them in two ranges separately, although the width is not long enough: one for the inertial-convective range 200 < r/η < 400 and the other for the range 30 < r/η < 60. Figure 29 compares the scaling exponents. ζ θ q in the inertial-convective range is larger than that of the viscous-convective precursor range. At the present level of resolution, no saturation of the exponents is found in the inertial-convective range, and the rate of increase of the scaling exponents in the viscous-convective precursor range becomes lesser and lesser as the order q increases and the asymptotic value 1.5 is approached [38,39,41,42,57,58,63,67]. The scaling exponents found in [45] are between the two cases in the present study. When we take global averages of the local scaling exponents over 30 < r/η < 400, the values are very close to those found by Chen and Cao [45], which suggests that their scalar scaling exponents in the inertial-convective range are underestimated.
We have seen that the mixed velocity-scalar structure function has a single exponent over the scaling range of 60 < r/η < 400 for fixed q. This is related to the constancy of the mean of the scalar flux in the global range. This motivates us to interpret ζ L q and ζ θL q as the scaling exponents of the energy transfer and the scalar variance transfer [48,56,59,65] i.e., |δu r | q = |r (r)| q/3 ∝ r ζ L q , |δu r δθ 2 r | q = |r θ (r)| q/3 ∝ r ζ θL q . ζ q q q/3 Chen & Cao (1997) Figure 29. Comparison of the scaling exponents for run 2. ICR is for the inertial-convective range 200 < r/η < 400, and VCPR is for the viscousconvective precursor range 30 < r/η < 60. C & C is the DNS data from Chen and Cao [45]. Figure 30 compares the scaling exponents ζ L q and ζ θL q for run 2 which were obtained by averaging the values over the range 100 < r/η < 300. Both curves pass through the point ζ L 3 = ζ θL 3 = 1; again we see ζ L q > ζ θL q . Also plotted are the curves for the She-Lévêque model with parameters γ = 2 3 , d 0 = 2 for the velocity; case 1: γ = 2 3 , d 0 = 1; case 2: γ = 2 3 , d 0 = 10 9 for the passive scalar. Agreement between the curves from the model for case 2 and from the DNS data is satisfactory, as reported by Lévêque et al [55,56]. The value d 0 = 10 9 close to unity suggests a sheet-like diffusive structure. The scaling exponents are listed in table 2.

Structure of the scalar and scalar dissipation fields
The energy and scalar dissipation rates (x) and χ(x) fluctuate in space and time. These fluctuations are considered to cause intermittency of the turbulence and scalar fields. It is therefore very interesting to see their spatial structures. Figures 31-33 show two-dimensional slices through the scalar field, the energy and scalar dissipation fields at the same time. The side of the square is 2π. The colour scale is determined by the following formula:  where A stands for (x, t) or χ(x, t). The motivation for using these transformations is to see the structure of the fields at both small and large amplitudes. If a linear transformation is used, high-amplitude events dominate the colour level, so that field patterns with weak amplitudes are difficult to observe. The structure in the passive scalar and its dissipation field are clearly seen. We strongly recommend that readers look at these figures by magnifying them on the online PDF version of this paper. There are large domains with the integral scale L θ in which the scalar amplitudes are   roughly constant. The boundaries between those regions are very sharp and highly convoluted. Also each large domain of high and low amplitudes contains smaller regions with weaker scalar fluctuations, and their boundaries are also sharp and convoluted. Corresponding to this, the scalar dissipation has large values at the boundary between the adjacent domains, as expected from the definition of the scalar dissipation. This can be well understood when one-dimensional profiles of the scalar amplitude and scalar dissipation are plotted in figures 34-36. In the figures, the scalar dissipation curve is shifted upward by a constant 5 for clarity. The scalar amplitude has a particular shape with a large scale plateau, sharp cliff and a deep valley, a mesa-canyon structure, which can be compared with the ramp-cliff structure in the scalar field in the presence of a mean scalar gradient. At the cliffs of the scalar field, large scalar dissipation occurs. At a first approximation, the set of points where the scalar dissipation is intense is complementary to the set of points where the scalar field is relatively smooth. The region of high scalar dissipation in figure 32 is long, continuous and wrinkled. The width and length of the high scalar dissipation domain are about 10η B and O(L θ ) respectively. Note that Sc = 1.
When (x) is compared with χ(x), the correlation between the two fields is weaker than that between the scalar and its dissipation. There are some domains in which the large structures in both fields resemble each other (upper left quarter), but also domains in which no resemblance is seen. This weak correlation can be seen also in small values of the correlation coefficients between the kinetic energy and scalar variance dissipations The latter value is comparable with the value C log log χ = 0.16 computed by Vedula et al [60]. The χ(x) field has a very thin structure compared with (x), and the intense domains of χ(x) are less space filling, indicating stronger intermittency of the passive scalar than the velocity. θ, χ x Figure 35. One-dimensional snapshots of θ(x, π, π) (red) and χ(x, π, π) (green) for run 2. The χ curve is shifted by 5 for clarity.

Resolution
To properly resolve the smallest scales of motion of the velocity and passive scalar fields, it is desirable for the value K max η = K max η B to be larger than unity. Sreenivasan [61] argued that since the intermittency produces a strong fluctuation of the energy and scalar dissipation, One-dimensional snapshots of θ(x, 3π/2, π) (red) and χ(x, 3π/2, π) (green) for run 2. The χ curve is shifted by 5 for clarity.
the smallest scales of the velocity and scalar are locally η min = (ν 3 / max ) 1/4 = η (¯ / max ) 1/4 , and η B,min = Sc −1/2 η min = η B (¯ / max ) 1/4 . Maximum value of the ratio max /¯ read from the PDF of by the present DNS is about 300 for run 2 (figure not shown). Thus we have η min = η B,min ≈ η /4, meaning K max η min = K max η B,min = 4 in the fairest resolution. And K max η = K max η B = 5 would be safer for a perfect resolution of the smallest scale because the present DNS tends to slightly underestimate max and χ max at the Reynolds numbers studied here.
There are two aspects about the spatial resolution of DNS. The first is that low resolution of small scales may affect the dynamics of the velocity and scalar fields, so that the statistics of the fields also are modified. To answer this question, one needs to compare the statistics computed by DNS using large K max η with those of low K max η, but this is a very expensive computation especially at the high Reynolds and Péclet numbers studied here. Our expectation is that the affected scales are those having several times η. However, the modification of the scalar dynamics due to K max η ≈ 1 can be inferred by comparing the present DNS results with those of previous studies. A comparison of the flatness of the scalar gradient with recent DNS and experiments, although some of them have a mean scalar gradient, shows that the scalar flatness is 14.0 for run 1 in the present DNS, 13.8 for Wang et al [46], 18.4 (perpendicular to the mean scalar gradient) for Yeung et al [20] and about 17 for Mydlarski and Warhaft [5,10] at R λ ≈ 200 and P λ ≈ 200. There is no inconsistency among the data up to the fourth-order moments, which strongly suggests that the scalar dynamics is not altered significantly. Of course, low-resolution effects may be found at higher orders, and more studies are necessary. Currently, there seem no systematic studies in this direction by using large and small K max η.
The second point is that when K max η and/or K max η B is equal to or smaller than unity, importance of the numerical error due to the low resolution depends on quantities and scales to be studied. For example, in the statistics of the scalar gradients, the second-order moments are affected only slightly, but the higher-order moments would be underestimated significantly. As found in figure 4 the dissipation spectra are contaminated at wavenumbers near kη ≈ 1, nevertheless the body of the dissipation spectra is well-resolved and statistically converges. The peak value 1.56 of the normalized scalar dissipation spectra is very close to 1.57 of that computed by Wang et al for which K max η = 2.96 and R λ = P λ = 195 [46]. This strongly suggests that when scales of interest are larger than a few multiples of η and/or η B , the low-resolution effects in the present DNS to the second-order statistics are negligible. When the order of the structure functions of the velocity and scalar increments becomes larger, the lower boundary of the range of scales over which the statistics is free from the numerical error increases. Although it is difficult to estimate the degree of error contamination due to the low resolution, the following arguments would be meaningful. Figure 37 shows the cumulative moments for run 2, where q = 10, z = δθ/ (δθ) 2 , P θ is the normalized PDF for the scalar increment. The symmetry of the PDF is assumed. The curves converge well for large z even at r/η = 2.8 = dx/η, and the PDFs with separation distance larger than the upper end of the diffusive range correspond to the curves with n 4. Figures 14 and 37 indicate that strong convergence of the structure functions for q = 10 and r/η > 10 is very encouraging for us to expect that the results in the present DNS are reliable as far as the statistics up to the 10th-order moments and at scales larger than about 10η are concerned.
Further quantitative studies about the resolution requirement of DNS are necessary. However, it is not possible currently for us to resolve both the inertial-convective and viscous-convective ranges simultaneously by using high-performance computers, not even using the Earth Simulator. One must choose the range to be analysed, and needs to compromise between the width of the inertial-convective range and the degree of faithfulness of the spatial resolution at scales near the finest grids. With this understanding, we have chosen to analyse the inertialconvective range in detail because detailed studies have not been performed in this range for DNS, and we have not quantitatively studied the statistics of the scalar derivatives.
To analyse properly the statistical quantities of turbulence, it is also required to obtain well-converged statistics in the sense of sample size, as seen in figure 4. Tails of PDFs are affected by sample size as well as by the spatial resolution. Usually, the time span of a few large eddy turnover times for time average is necessary to obtain smooth curves of PDFs and structure functions at higher order. This is relatively easier to do when compared to satisfying the spatial resolution requirement. The length of time span of runs 1 and 2 exceeds two large eddy turnover times.

Summary and discussion
We have numerically studied the statistics of a passive scalar convected by incompressible turbulence when Sc = 1. The passive scalar spectrum in the inertial-convective range was found to be E θ (k) = C OC¯ −1/3χ k −5/3 over a finite range of wavenumbers and the Obukhov-Corrsin constant C OC was found to be C OC = 0.68 which is consistent with experimental values. The spectral bump is more significant for the scalar than the velocity.
Yaglom's 4 3 -law for the mixed velocity-scalar increment correlations is approached when the Péclet number is increased. It was argued that the 4 3 -law holds over a wider range of scales than the 4 5 -law for the velocity triple correlation when Sc 1. This is because the average scalar transfer flux has the same constant value throughout the inertial-convective and the viscousconvective ranges. The structure functions for the velocity and passive scalar increments were computed and compared. It was found that the structure functions for the passive scalar have two different scaling ranges, one for the inertial-convective range and the other for the viscousconvective precursor range. Although the Schmidt number 1 is not large, the non-locality of the scalar transfer in the scale space is efficient for a narrow viscous-convective precursor range to be observed.
The local scaling exponents of the velocity, scalar and mixed velocity-scalar structure functions were computed as functions of the separation distance. In the inertial-convective range, the scaling exponents of the passive scalar structure functions are found to be increasing functions of the order and smaller than those of the velocity. But it is not certain whether saturation of the exponents occurs or not when R λ and q increase. On the other hand, those in the viscousconvective precursor range appear to be saturated at about 1.5 as the order increases.
The structure function of the mixed velocity-scalar correlation, which is interpreted as the flux of the passive scalar multiplied by r, behaves smoothly as if there was a single scaling range of the separation r over the inertial-convective range and viscous-convective precursor range. This is a new finding of the present high-resolution DNS. The local scaling exponent maintains a constant value for each order, and is smaller than that of the energy transfer flux for q > 3, indicating stronger intermittency of the scalar transfer flux.
More accurate determination of the scalar scaling exponents in the inertial-convective range requires a longer scaling range, or equivalently higher Péclet numbers; we need larger DNS for the passive scalar. However, it is certain that the intermittency of the passive scalar is stronger than that of the velocity and the scalar dissipation has a more singular structure than the energy dissipation. The mesa-canyon (or ramp-cliff) structure is generic for the passive scalar convected by the turbulence, regardless of whether the turbulent velocity field is synthetic or generic, and whether or not a mean scalar gradient exists. These structures suggest a persistent interaction between macro-and microscales of motion, which can be identified as a key ingredient in understanding the intermittency of the scalar.
When compared with the velocity statistics, the fluctuation level in the passive scalar is stronger, which reflects the stronger intermittency of the passive scalar. The interactions in wavenumber space are more non-local than the interactions in the velocity field. Even when Sc = 1, this non-locality manifests itself as a precursor to the viscous-convective range. Since there is no pressure term in the equation governing the passive scalar, the scalar field suffers from a large scale deformation by convective action alone, and the local isotropization and mixing by the pressure which are so important for the velocity field are entirely absent. Thus, the scalar can be expelled from the body of the vortical motion and piled up at the periphery or at stagnation points, leading to the mesa-canyon structure. These facts suggest that scalar statistics converge more slowly than velocity statistics, requiring larger sample sizes and wider scale separation to obtain accurate and reliable data.
It is also interesting to compare the passive scalar turbulence with the Burgers turbulence. Although there is a difference between them in that the velocity field of the Burgers turbulence is not passive, they share common features that both have no pressure term, ramp-cliff (mesacanyon) structure, and localized structure (thin shell structure in three dimensions) of high dissipations [62]. As an outcome of these similarities, the intermittency of two fields is stronger than that of the incompressible fluid turbulence. At the same time, strong intermittency puts more severe conditions on the spatial resolution and statistical convergence required for their DNS studies. Therefore, DNS study of the passive scalar is more difficult than in the case of the velocity, and remains a major challenge to turbulence theory and computational science.