Turbulent mixing: the roots of intermittency

Mixing in fully developed incompressible turbulent flows is known to lead to a cascade of discontinuity fronts of passive scalar fields; these fronts make the tracer fields strongly intermittent and give rise to an anomalous scaling of structure functions. Here, a simple model, essentially a one-dimensional (1D) variant of Baker's map, is developed, which captures the main mechanism responsible for the emergence of the discontinuities. The model is studied both numerically and analytically. In particular, the structure function scaling exponent ζp is derived; for the Kolmogorov turbulence, . The analytical findings are consistent with simulations, and explain the results of a series of numerical and experimental studies.


Introduction
Fully developed turbulence is known to be highly intermittent, characterized by non-Gaussian statistics and non-vanishing probabilities of extreme events. The basic equations defining the evolution of the turbulent media are relatively simple; therefore, the turbulence is a convenient polygon for acquiring knowledge of intermittency. For instance, the concept of structure functions, introduced for describing the non-Gaussianity of turbulent fluctuations [1], has evolved into a nowadays very popular multifractal formalism. So, progress in understanding the driving forces of intermittency is of great significance beyond turbulence theory itself. In this regard, the best research polygon is provided by the turbulence of passive scalars. Indeed, as compared e.g. with hydrodynamic turbulence, the turbulence of passive scalars is characterized by stronger intermittency and simpler evolution equations.
The empirical aspects of the intermittency of passive scalar turbulence are rather well understood; also, there has been great progress in theory (cf [2]- [6]). Regarding the structure function scaling exponents ζ p , there are detailed experimental and numerical studies (cf [7]- [13]). In particular, it has been shown that the anomalous scaling is the result of the formation of tracer density discontinuity fronts [12]. There are also hierarchical models [14]- [16] and models based on the assumption of the multifractal nature of the passive scalar field [17,18], which can fit the experimental data relatively well. However, they rely on a couple of phenomenological hypotheses, and therefore their usefulness for understanding the origins of intermittency is limited. As for analytic studies, several important rigorous results have been obtained, cf [2,19,20]; still, the explanation of the experimental curves of ζ p is not complete.
Here, we develop a novel, simple model for the turbulent mixing of tracers, which is in essence a stochastic hierarchical Baker's map model. Baker's map performs first an area-preserving affine transform on a region, and then 'folds' it, to restore the initial shape of the region. It mimics the stretching-folding action of chaotic fluid motion, and is therefore well suited for the modelling of chaotic mixing, cf [21]. We modify the mapping so that it would not introduce discontinuities; the mapping is applied hierarchically to different-sized subregions of the polygon.
Firstly, we provide a physical motivation of the model. Further, we analyse theoretically the properties of our model, and present the simulation results. The theoretical results are in good agreement with our numerics, as well as with the earlier numerical and experimental studies of the ζ p -curves. Therefore, we argue that our model captures the main mechanism responsible for the intermittency and anomalous scaling of the passive scalar turbulence. Still, it should be emphasized that the model is too robust for explaining all the features of the ζ p -curves (such as e.g. the differences between two-dimensional (2D) and 3D flows).

Motivation of the model
In what follows, we assume that the tracer density θ(r, t) evolution is described by a simple diffusion equation, where v(r, t) is a turbulent velocity field in a 3D or 2D space, f (r, t) is a forcing, and the seed diffusivity κ is assumed to be very small (but not zero). Finally, it is assumed that the Characteristic time of eddies of size a scales as τ ≈ a 2−ξ . Due to the combined effect of large and small eddies, low-and high-density regions (black and white, respectively) are brought into contact within a finite time.
velocity field is statistically stationary, homogeneous and scale-invariant; we address both timecorrelated flows and flows with Kraichnan statistics (delta-correlated in time).
It should be noted that we neglect the intermittency of the underlying velocity field. According to the experimental and numerical evidence, the passive scalar intermittency is stronger than the velocity field intermittency (e.g. characterized by greater anomaly 2ζ 2 − ζ 4 ), cf [4]. Thus, one can expect that the former dominates over the latter, and the behaviour of tracers in real (intermittent) velocity fields is very similar to that in idealized Gaussian fields (for numerical evidence, cf [22]).
The mixing effect of turbulent flows is most intuitively characterized by the growth of the distance between two tracer particles r . In the case of Kraichnan statistics, where ξ is the velocity field smoothness exponent, defined as the scaling exponent of the nonconstant part d i j of the velocity field correlation function v i (r, t)v j (r , t ) . So, with a proper time unit, the distance doubling time is estimated as τ ≈ r 2−ξ ; cf [2]. Equation (2) holds also for a wide class of time-correlated flows, assuming that an appropriate value of ξ is used (for a more detailed discussion of this matching, cf [23]); in particular, the Kolmogorov flows are characterized by τ ≈ r 2/3 and therefore are matched with ξ = 4/3. For the sake of simplicity, we consider 2D geometry; generalization to the 3D geometry is straightforward. The formation of tracer discontinuities can be qualitatively explained as follows. Firstly, we decompose the velocity field into components of different-characteristics space scale a, v a (r, t) = a |k| −1 <2a v(k, t) e ikr dk, where v(k, t) is the Fourier transform of v(r, t). The characteristic timescale τ a is defined as the time needed for the field v a (r, t) to transport a tracer particle to a distance of the order of a; according to equation (2), τ a ≈ a 2−ξ .
Suppose that initially (t = 0), there is a constant tracer gradient: θ(r, 0) ≡ x. In figure 1, the regions θ 0, 0 < θ 1 and 1 < θ are shown in black, grey and white, respectively. The minimal distance between the isolines θ = 0 and θ = 1 will start decreasing (they are turbulently stretched, and the surface area between them is conserved). Initially, this approaching is dominated by the largest eddies fitting between these lines, i.e. the eddies of approximately unit diameter, characterized by τ 1 ≈ 1. After a unit time, the distance between some segments of the isolines will get reduced approximately by a factor of two. From now on, the distancedecreasing rate will be dominated by two times smaller eddies, and the characteristic timescale is reduced by a factor of 2 2−ξ . The process will continue ad infinitum, leading to the contact of segments within a finite time (the timescales form a geometric progression). This scenario holds both for real flows and for Kraichnan flows. In fact, there are differences in the character of the isoline transport. In the case of realistic time-correlated flows, the motion of the isolines, due to the field v a (r, t), is convective; the characteristic time τ a is estimated as the ratio of the scale a and the characteristic velocity v a . For the Kraichnan flows, the motion of the isolines is of a diffusive character; the characteristic time τ a is estimated as the ratio of a 2 and the diffusivity d i j ≈ a ξ . However, the end result (i.e. the shape of an isoline upon passing a time τ a ) is qualitatively the same in both cases (sketched in figure 1).
We aim to construct a model that mimics the evolution of the tracer density profile along the x-axis (any 1D cross section). To begin with, we consider only the effect of an 'a-flow' v a (r, t). In incompressible flows, the evolution of scalar fields is driven by the stretching-folding motion of the fluid (cf [24]). Such a stretching-folding motion is provided by a simple shear flow, as depicted in figure 2. Now, consider the tracer density profile along the x-axis in figure 2: the initially monotonous curve θ (x) is replaced by a new profile θ (x) with a 'kink'. The kink emerges because a descending segment is replaced by the sequence of descending, ascending and descending segments. In the idealized version, all these curve segments are mirror images of each other, and are the result of a 3-fold 'compression' of the initial curve segment along the x-axis. Such a mapping M a,c : θ (x) → θ (x), is essentially a variant of Baker's map, figure 3(A); here, a denotes the size of the vortex and c, its midpoint.

The model
Mapping M a,c models the basic elementary process responsible for the emergence of the discontinuity fronts: it mimics the effect of a single shear flow on the tracer density profile. As can be seen from figure 1, for a fully developed turbulence, the fronts emerge as a result of the collective action of different-sized vortices. So, the mixing action of the turbulent fluid is modelled by an iterative application of the mapping M a t ,c t to some initial profile θ 0 (x) with random values of the parameters a t and c t : note that, here, t plays the role of discrete time.
In order to match statistically homogeneous turbulence, the probability density function (PDF) of this mapping over the parameter c needs to be homogeneous. Also, our model needs to match the stretching statistics of the velocity field (2). To this end, the PDF of the mapping over the parameter a needs to be such a power law that the waiting time T a between two subsequent mappings of size a at a site x = c (i.e. satisfying the conditions a t ∈ [a, 2a] and c ∈ [c t − a t /2, c t + a t /2]) would scale as T a ≈ a 2−ξ . Now, let us turn to the question of how to model the boundary and initial conditions. Initial conditions in the form θ t=0 = θ 0 (r) can be modelled in a straightforward way, by the initial profile θ 0 (x). The easiest way to model the situation when there is no tracer flux through the container boundaries is by applying the periodic boundary condition θ t (x + 1) ≡ θ(x). The effect of forcing in equation (1) can be mimicked also in a very straightforward way, by additional iterations Finally, let us consider the boundary condition θ x=0 = 1. If θ is the temperature, this corresponds to a heating wall with a fixed temperature. Notice that if a large vortex is positioned close to the 'hot' wall, it provides a persistent 'washing' of the wall, and hence, a persistent heating of the surface layer of the fluid to the wall temperature. The heating effect (i.e. the volume of the heated fluid) is proportional to the volume of the vortex. Such a process can be matched as depicted in figure 3(B). Assuming that the centres of the vortices are evenly distributed over the whole volume (c ∈ [0, 1]), the edge of the mapping can extend beyond the boundary of the fluid. If this is the case, for instance if c t < a t /2, the profile θ t (x) is extended to the region x < 0 with the boundary value. Then, the mapping is applied in the same way as described earlier (with a single modification: the 'compression' factor is increased from 3 to 3( 3 2 − c t /a t ), to fit the entire 'vortex' inside the region x > 0; see figure 3(B)).
To conclude, we need to discuss how to model the effect of the seed diffusivity. For any nonzero diffusivity κ, the diffusion smoothes the tracer density fluctuations at a microscale δ, for which the effective Peclet number P δ ≈ 1. From the equality of diffusion and mixing times, τ δ ≈ δ 2−ξ ≈ δ 2 /κ, we obtain δ ≈ κ 1/ξ . In order to take into account such a smoothing, the mapping M a t ,c t is modified so that, apart from the effect depicted in figure 3, it includes also averaging over a sliding window of width δ.
For numerical simulations, averaging over a sliding window can be avoided if the dissipation scale δ is used as the discretization scale. Then, the tracer density profile is stored as an array θ i,t ≡ θ t (iδ), i = 1, 2, . . . , N with N = κ −1/ξ . Now, if the mapping transforms a point i into a point j, the iteration formula provides an effective averaging over the dissipation scale.

Theoretical analysis of the model
Our scaling analysis is based on the PDF f a ( a ) of the difference a (x) = |θ a (x + a 2 ) − θ a (x − a 2 )| between the mean values of the tracer densities for neighbouring segments of length a. Here, the local average is defined asθ a (x) = 1 a x+ a 2 x− a 2 θ(y) dy. To begin with, we consider what will happen if segments AB and BC, characterized by mean densitiesθ AB andθ BC , are transformed by a mapping, see figure 3(C). The segment AB is transformed into three thrice smaller segments, two of which are marked as D E and G H ; the mean densities for these segments are equal to that of the segment AB. The same applies to the segments BC, E F and F G. So, the density drop between the segments D E and E F is equal to that between the segments AB and BC, i.e. to a (B). However, there is no density drop between the segments E F and F G. Apparently, the density drop a/3 (x) takes all the intermediate values between a (B) and 0, as x moves from E to F; the dependence on x is approximately linear (at least in the neighbourhood of F). Consequently, as a result of the mapping, the probability f a ( a ) d a , associated with the density drop a , contributes to the PDF f a/3 evenly over the range of values a/3 ∈ [0, a ]. This mechanism relates all the values to the values of f a ( a ) (because mappings are continuously being applied to the profile θ (x)); mathematically, (assuming that the maximal value of is 1). Note that a similar integral equation has been used to describe the mixing of a tracer emitted from a point source [25]. It should be emphasized that equation (3) is obtained using two implicit assumptions.
(i) We do not consider the effect of those mappings, the size of which is either significantly larger or smaller than a and b. It can be argued that the effect of smaller-size mappings  (3)). So, the probability f a ( ) d associated with the drop at scale a (denoted by the leftmost black box) is redistributed over a range of smaller values of at thrice smaller scale; this is symbolically marked by arrows at step 1 * of the scheme. Further, if any of the now thrice smaller segments is hit by a much larger mapping (i.e. with a a 1 /3), the existing density drop is just reproduced at thrice smaller scales (the resulting flow of the PDF is depicted by horizontal arrows at step 2). If we want to relate the PDF at scale a 2 to the PDF at scale a 1 , we need to count the number of convolutions n, i.e. the number of steps marked by ' * ' in this scheme. Then, it is convenient to introduce the effective compression factor k, defined as the average distance in logarithmic scale between subsequent convolutions (so that ln k = 1 n ln(a 1 /a 2 )).
is insignificant at our scale, because they preserve the average densityθ b . This is true, if the mapping falls entirely into the segment; if it falls at the edge,θ b will be changed, but the change remains relatively small. However, for very small values of ξ < ξ 0 , when small vortices are much more frequent than the large ones, this argument will no longer be valid. (ii) We can neglect the effect of larger vortices. This is actually not true: larger-size mappings compress the profile without reducing the density drop via the process depicted in figure 3. Such a process corresponds to a direct transfer f a ( ) → f a/3 ( ), without the convolution in equation (3). So, on average, the profile will be compressed more than thrice before entering the convolution stage. Hence, the effect of larger vortices can be taken into account by using an effective, somewhat increased compression factor k = a/b > 3. The concept of the effective compression factor is illustrated in figure 4.
Recall that we expected k 3; this inequality is not satisfied for ξ < 1. So, for ξ < 1, the assumption (i) is not satisfied, i.e. ξ 0 = 1. The result ξ 0 = 1 is directly applicable only to our 1D model, when the compression factor is always k = 3; in the case of real 2D or 3D turbulence, k may take different effective values and, hence, the critical value ξ 0 may depart from 1.

Numerical analysis
We have implemented the above described model (with boundary conditions θ x=0 ≡ 1 and θ x=1 ≡ 0) numerically for several values of ξ . The array length (which corresponds to the ratio of the tank width and diffusion scale) was taken to be equal to N = 10 6 . For each value of ξ , the observation time was long enough to include at least 200 full decorrelations of the field θ t,i . Simulation results are presented in figures 5 and 6. The small mismatch between the curves and data points in the left-hand side of figure 5 can be explained by finite-size effects and somewhat poorer statistics of extreme events. In the case of figure 6, one should recall that our theoretical analysis (and, hence, curve (a)) is applicable only to the range ξ 1.

Discussion of the results
It should be emphasized that our model is very simplified. In fact, it is difficult to expect that a 1D model can capture all the details of the 2D or 3D reality. So, there is numerical evidence that3Dflowsarecharacterizedbyaslightlylargerscalinganomaly2ζ 2 − ζ 4 than 2D flows, see figure6 [30].Meanwhile,westatedthatthebasicmechanismresponsiblefortheintermittencyis (b) the Kraichnan formula [28]; (c) perturbation theory [29]; (d), (e) the numerical results of 2D and 3D Lagrangian simulations [19,30]; (f) simulations with our 1D model. the same both for 2D and 3D geometries, and both will be modelled by the same hierarchical 1D Baker's map. Now, comparing the predictions of our model with the experimental and numerical data (figures 5 and 6), we can conclude that the initial claim remains valid: both 2D and 3D data are matched fairly well. However, the model is not detailed enough to explain the empirically observed small difference between the 2D and 3D flows.
Perhaps a more important issue (at least from the theoretical point of view) concerns the asymptotic behaviour of the structure function scaling exponents ζ p at the limit p → ∞. It has been believed that there is a saturation of these exponents, cf [2], which is in contradiction to the logarithmic growth predicted by equation (4). This belief is based on the following facts. Firstly, it has been shown that such a saturation indeed exists at the limit of large dimensionalities d of the geometrical space, for d(2 − ξ ) 1 [20]. Secondly, some experiments seem to support the saturation, cf [12]. In particular, it has been pointed out that if there is a saturation, then there are also so-called mature fronts-the fronts across which the tracer density drops by a non-small fraction (e.g. 50%) of the largest density difference for the whole polygon; the existence of the mature fronts has been confirmed numerically [12].
When discussing the issue of the logarithmic growth versus saturation, it should be stressed that the logarithmic growth is very slow: in practice, the difference is so small that e.g. it is difficult to make conclusive judgments based on the experimental data, see figure 5 (one should also bear in mind that the uncertainties of ζ p grow rapidly for larger values of p). So, we can draw two conclusions. Firstly, independently of the question of which asymptotic behaviour appears to be the correct one, our model remains a reasonable approximation of the real tracer dynamics. Secondly, the existing numerical and experimental data are not precise enough to exclude one or another asymptotic form.
Finally, we need to comment on the theoretical results for d(2 − ξ ) 1. As mentioned above, the basic mechanism responsible for the tracer intermittency is independent of the dimensionality d 2. However, higher dimensionalities are characterized by slightly stronger intermittency (i.e. by smaller values of ζ p for large indices p). So, the combination of two limits, d 1 and p 1, is a special case, which was not intended to be described by our model. Hence, there are two possibilities: firstly, our model is not detailed enough to describe correctly the limit p → ∞ ; secondly, it does describe correctly this limit, but only for moderate dimensionalities d. Further studies are needed to resolve the question of which option is valid.

Conclusion
Although simplified, our model captures the basic mechanism responsible for the intermittency and anomalous scaling of the passive scalar turbulence. It explains, with a reasonable accuracy, the results of previous experiments and simulations. The model also opens perspectives for a wide range of further studies. Firstly, there are some open questions, e.g. the incompatibility of equation (4) with [20]. Answers to these questions may be obtained by adding more details to our model and analysing the robustness of our results. Perhaps the simplest addition would be to decrease the symmetry of the mapping in figure 3(A) by making the lengths of the three segments (forward, reverse and forward ones) unequal. Secondly, the model presented above can be used as an efficient tool for studying various problems related to the turbulent transport and mixing, both numerically and analytically. In particular, we anticipate further progress in understanding the rain initiation process in warm clouds [31]. An important stage of the rain initiation process is to widen the droplet size spectrum by mixing the air volumes of different droplet sizes. We also expect that the model can help in understanding the mechanisms creating the small-scale anisotropy of the passive scalar fields [4]. Further, it is tempting to study the mixing of passive scalars in the case of point sources, cf [25], as well as to extend our approach to the modelling of kinematic dynamos, cf [32].