Eddy diﬀusion
December 7, 2000

1

Introduction

Perhaps you have heard that turbulence is the most diﬃcult problem in
ﬂuid mechanics and, according to some, the greatest unsolved problem inphysics. One indication of the diﬃculty is that it is impossible to give asatisfactory deﬁnition of a “turbulent ﬂow”. But everyone agrees that one
property of turbulence is greatly enhanced transport of passive contaminants.
For example, relying only on molecular agitation, a dissolved sugar moleculetakes years to diﬀuse across a coﬀee cup, and on that time-scale the coﬀee
will surely evaporate. With a spoon the coﬀee drinker can create eddies thattransport dissolved sugar throughout the cup in less than one second. Thisis an example of eddy diﬀusivity.
Fluid mechanics textbooks often often justify eddy diﬀusivity by appealing to an analogy between turbulent eddies and molecular diﬀusion — perhaps this notion originates with G.I. Taylor’s 1915 paper entitled “Eddymotion in the atmosphere”? In any event, the molecular analogy, supplemented with some hand-waving, leads to the notion of an eddy diﬀusivity
and for many scientists this is the end of the turbulence problem.
Our goal in this lecture is to explain very explicitly the assumptions
behind Taylor’s “proof by analogy” and to illustrate the interesting points at
which the analogy fails. We will pursue this program by working with somevery simple model ﬂows for which analytic results, such as expressions forthe eddy diﬀusivity, are available. As you will soon see, these model ﬂows
do not greatly resemble turbulence, but then neither does molecular motion!
Our excuse is that soluble examples are always diverting and educational.

1

2

The renovating wave model

A recipe for constructing soluble models
The main problem in analyzing transport is solving the diﬀerential equations
which describe the motion of particles in even very simple ﬂows. However
there is a class of ﬂows for which this task is trivial. These are steady andunidirectional ﬂows, such as u = sin y. A particle which starts at (a, b) at
t = 0 ﬁnds itself at (a + τ sin b, b) at t = τ . This is dull, but it becomes moreinteresting if at intervals of τ we “renovate” the ﬂow by randomly picking a
new direction along which the velocity acts. In this way we can constructa sequence of iterated random maps and calculate diﬀusivities, and otherstatistical properties, by averaging the exact solution. I learned of this trickfrom the literature on dynamo theory. The book Stretch, Twist, Fold: theFast Dynamo is highly recommended for students interested in all aspects ofstirring and mixing. .

2.1

The renovating wave (RW) model

As a particular example we now formulate the renovating wave (RW) model.
We divide the time axis into intervals
In ≡ {t : (n − 1)τ < t < nτ } ,(1)and in each interval we apply a velocity, u = (−ψy , ψx ), derived from thestreamfunction
ψn (x, y, t) = k −1 U cos[k cos θn x + k sin θn y + ϕn ] ,(2)where θn and ϕn are independent random variables uniformly distributed in
the interval [−π, π]. Thus in each In there is a steady, unidirectional velocity
with sinusoidal proﬁle (a single wave). There is sudden and complete loss of
all information about the past velocity at t = nτ because at these instants
we “renovate” the velocity by picking new random angles θ and ϕ. (Thismeans that the velocity correlation function, C(t), is zero if t > τ .)The renovating wave model can be nondimensionalized by using k −1 as aunit of length and and 1/(U k) as a unit of time. With this choice, the model
contains a single dimensionless parameter, τ∗ ≡ τ kU , which is a measure of2the persistence of the motion. Much of the literature on random advectiondiﬀusion uses model velocity ﬁelds which are δ-correlated in time. We canrecover this limit as a special case by taking τ∗ → 0.
Using dimensionless variables, a particle which is at xn = (xn , yn ) at
tn = nτ∗ moves to xn+1 at t = (n + 1)τ∗ where(xn+1 , yn+1 ) = (xn , yn ) + τ∗ sin (cn x + sn y + ϕn ) (sn , −cn ) .

(3)with sn ≡ sin θn and cn ≡ cos θn . Thus motion in the renovating wave
problem is equivalent to an iterated sequence of random maps.

2.2

The single–particle diﬀusivity

It is very easy to calculate the diﬀusivity in the RW model (and much morediﬃcult to interpret the answer). The average of a function of the two random
angles θ and ϕ (suppress the subscript n) is deﬁned by


dϕ
dθ
f  =f (θ, ϕ).
(4)2π
2π
Therefore, using (3),
(xn+1 − xn )2  =

τ∗2
.
4

(5)

The computation is trivial if the integral over ϕ is evaluated ﬁrst.
In (5), following our previous discussion based on Einstein’s derivationof the diﬀusion equation, we are computing the statistics of dispersion along
the x-axis. Because the renovating wave model is isotropic, dispersion in they-direction is identical to that in the x-direction.
Because all of the waves are independent and identically distributed itfollows that after n renovation cycles
τ2
(xn − x0 )2  = n ∗ .
4(6)

But t = nτ∗ , and (xn − x0 )2  = 2Dt, so that using dimensionless variables
the diﬀusivity is
D=
3

τ∗
.
8

(7)

t=1τ

t=2τ

t=3τ

t=4τ

t=5τ

t=6τ

t=7τ

t=8τ

t=9τ

t=10τ

t=11τ

t=12τ

Figure 1: Stretching of a small spot, r
1 where r is the initial radius of thespot, by a succession of random sinusoidal ﬂows. The dotted circle is the initial
spot.
Sometimes D is referred to as the single-particle diﬀusivity. “Single-particle”
emphasizes that D strictly applies only to the RMS displacement of a particle
from its initial position; D contains no information concerning the deformation of a patch of tracer, nor of any other quantity involving correlated motion. Thus, using dimensional variables, the diﬀusivity in (7) is D = U 2 τ /8,
which is independent of k. Because D is independent of the scale of the wave,
even a spatially uniform, but random-in-time velocity (the case k = 0), hasa single-particle diﬀusivity.

4

t=1τ

t=2τ

t=3τ

Figure 2: Stretching of a blob with r = 1, where r is the initial radius. The dotted
circle is the initial patch.

2.3

Deformation of variously sized blobs

To emphasize the importance of understanding more than single-particle diffusivities we take a digression and illustrate how the deformation of an initially circular blob of ﬂuid depends on the blob radius r. (Recall that wehave used k −1 as unit of length; in terms of dimensional variables the relevant
nondimensional parameter is kr.)
If the initial blob is much smaller than the wavelength of the velocity
then on the scale of the blob the velocity proﬁle is a linear function of thecoordinates. Because of this simplicity, the ﬁrst few iterations deform thecircular blob into an ellipse which must have the same area as the initialcircle. We will see in the next lecture that the major axis of the ellipse grows
exponentially while the minor axis shrinks so that the area is ﬁxed. Oncethe dimensions of the ellipse are comparable to the wavenumber of the ﬂow,
more complicated deformations occur. Ultimately the blob will be stretched
into a folded ﬁlament as in ﬁgure 1.
The blob has the same scale as the velocity ﬁeld if r ∼ 1. Because there
is no scale separation there is no easy description of the action of the ﬂow on
the blob, see ﬁgure 2.
If r
1 then we are in the “eddy diﬀusivity” limit in which the scaleof the velocity ﬁeld is much smaller than the scale of the tracer. This caseis shown in ﬁgure 3. The action of the waves perturbs the edge of the blob,
making it look “fuzzy”. In fact, the area is preserved, but the circumferenceof the blob grows exponentially. We will be discussing this type of problem
for the remainder of this lecture.

5

t=1τ

t=2τ

t=3τ

t=4τ

t=5τ

t=6τ

t=7τ

t=8τ

t=9τ

t=10τ

t=11τ

t=12τ

Figure 3: Stretching of a big blob r

1, where r is the initial radius of the blob.
The dotted circle representing the initial patch may not be visible beneath the
wiggly boundary of the blob.

2.4

The Lagrangian correlation function

In (7) we gave the diﬀusivity of particles moving in an ensemble of renovatingwaves. How do we obtain the Lagrangian velocity autocorrelation function
and verify Taylor’s formula that
 ∞D=C(t) dt ?
(8)06

Considering this question, we encounter an annoying technical diﬃculty: ourderivation of (8) assumes that the velocity statistics are stationary. But therenovating wave ensemble, as we deﬁned it back in (1) and (2), is not a
stationary stochastic process. This is because with our original deﬁnition
all members of the ensemble renovate at the same instants t = τ , t = 2τ
etcetera. In order to obtain a stationary process we should initiate diﬀerent realizations at uniformly distributed points during the renovation cycle.
Thus, for realization number j, we pick a random time τ (j) which is uniformly
distributed in the interval [0, τ ] and renovate ﬁrst at t = τ (j) and then subsequently at t = τ (j) + τ , t = τ (j) + 2τ etcetera. With this new and improved
formulation of the RW model the Lagrangian correlation function of u(t) isa “triangular” function:U2tC(t) =1−
H(τ − t) ,(9)4
τwhere H is the step function and U is the velocity in (2). The area of under
this correlation function is D = U 2 τ /8.

3

The eddy diﬀusion equation

3.1

The ensemble averaged Green’s function

Now that we have obtained the RW diﬀusivity in (7) we turn to the derivationof the eddy diﬀusion equation. For each realization we introduce the Green’s
function which is
Gt + u·∇G = 0,

with G(x, x0 , 0) = δ(x − x0 ).

(10)

The solution of the problem above isG(x, x0 , t) = δ(xt − x0 ) ,(11)where xt is the position at time t (in a particular realization of u) of theparticle which started at x0 .
The ensemble averaged Green’s function isg(r, t) = G(x, x0 , t) ,7r ≡ |x − x0 | ,(12)where we have assumed that the random velocity is isotropic, homogeneousand stationary so that g can depend only on the distance r and the elapsed
time t.
Possessing g(r, t), we can then represent the ensemble-averaged solutionof the initial value problem
ct + u·∇c = 0,
as the convolutionc(x, t) =c(x, 0) = c0 (x) ,(13)c0 (x − x )g(|x |, t) dx .

(14)(We are assuming that the initial condition c0 is the same for all realizations.)
At this point, the analogy between (14) and the master equation of lecture
1 is obvious. With the master equation in mind, we can anticipate that avariant of Einstein’s derivation of the diﬀusion equation can be applied to(14). Rather than develop a general derivation we prefer to use the renovating
wave model as a concrete illustration of how one can obtain g, and then passfrom the integral equation in (14) to an approximate diﬀusion equation.

3.2

The averaged Green’s function of the RW model

There are at least two ways of obtaining g(r) in (12) for the RW model: thehard, straightforward way (see the appendix) and the easy, devious way. Let
us be devious.
We begin by calculating the probability density function (PDF) of displacements in a single pulse of the RW model. Because the ensemble ofvelocities is isotropic and homogeneous there is no harm in supposing thatthe particle is at the origin and the x-axis is aligned with the direction of
the velocity. That is, put (xn , yn ) = (0, 0) and θn = π/2 in (3). Thus, thedisplacement r produced by a single pulse isxn+1 − xn = τ∗ sin ϕn ,and r = |xn+1 − xn |.

(15)

The PDF of the random variable r can be obtained from the PDF of ϕ, thatis P (ϕ) = 1/2π, using the rule for transforming probabilities:
 dϕ 

2 H(τ∗ − r) 
P (r) =
P (ϕ)  , ⇒ P (r) = 
.
(16)
dr
π τ∗2 − r2
8

4000

3000

2000

1000

0

0

0.2

0.4

0.6

0.8

1

r

Figure 4: Displacements of 40, 000 particles in independent realizations of thethe RW model. The left panel shows the ﬁnal position of particles which all startat the center of the the circle after one pulse of the wave. The density of points
corresponds to g(r) in (17). The histogram on the right shows the number ofparticles at a distance r from the center; this is the function P (r) in (16).

In (16) H(τ∗ − r) is a Heaviside step function which ensures that there are
no displacements greater than τ∗ . (The sum in (16) is because there are four
values of ϕ corresponding to a single value of r.)
The averaged Green’s function is now given byg(r) =

P (r)
,
2πr

⇒

g(r) =

1 H(τ∗ − r).
π 2 r τ∗2 − r2

(17)

The geometric factor 2πr is included because g(r) is a concentration. That
is, P (r)dr the expected number of particles which fall into the diﬀerential
annulus between r and r + dr and g(r) is the expected number of particles
per area in this same annulus; see ﬁgure 4.
Now that we have the averaged Green’s function of a single pulse we canobtain the evolution the ensemble averaged concentration, c, over many
pulses. Because each pulse is independent of the preceeding pulses we havec (x, (n + 1) τ∗ ) = c(x − x , nτ∗ )g(|x |) dx .
(18)
The master equation above, with g(r) in (17), is an exact description of theevolution of c under advection by the RW model.
93.3

The diﬀusion equation

With the master equation (18) in hand, we can use Einstein’s approximations
to obtain the diﬀusion equation. Using the dimensionless variables of therenovating wave model, we haveτ∗ct ≈ ∇2 c .
(19)
8
We leave this as a homework exercise and instead we take a diﬀerent routeto (19).
Because the Fourier transform of a convolution is the product of theFourier transforms, we can simplify (18) by transforming. The Fourier transform of f (x) is deﬁned here1 as


1
−ik·x
˜
(20)
f (k) = ef (x) dx ,f (x) =
eik·x f˜(k) dk .
2π
Applying the transform to (18) we obtainc(k, nτ∗ ) = g̃(k)n c̃0 (k) ,k ≡ |k| .

(21)With a good table of integrals one can discover that the Fourier transformof the averaged Green’s function, g(r) in (17), isg̃(k) = J02 (kτ∗ /2) ,(22)

where J0 is the Bessel function.
The diﬀusion equation describes the evolution of large spatial scales,
which is the same as small wavenumbers. This means that we simplify (21)by taking kτ∗ /2
1 and using the approximation J0 (kτ∗ /2) ≈ 1 − (k 2 τ∗2 /16)
to writec(k, nτ∗ ) ≈ exp n ln 1 − (k 2 τ∗2 /8)c̃0 (k) .

(23)

But now, since n = t/τ∗ and ln[1 − (k 2 τ∗2 /8)] ≈ −k 2 τ∗2 /8, we havec(k, t) = e−Dk t c̃0 (k)2(24)

where, as in (19), D = τ∗ /8. Equation (24) is the equivalent to the decay of
Fourier components given by (19).
This derivation based on Fourier analysis explicitly recognizes that thediﬀusion approximation is valid only for wavenumbers which satisfy kτ∗ /2
1. This is a precise statement of the scale separation assumption which
underlies Einstein’s approach.
1

By denoting the wavenumber with k we are recycing notation used in (2).

10

4Ensemble averages and single realizations

In hydrodynamic dispersion, particles which begin at neighbouring points
have similar histories in any single realization. Marbled endpapers in oldbooks were produced by ﬂoating coloured inks on water, stirring the surface, and then capturing the swirls by carefully lowering a sheet of paper
onto the inky ﬁlm [3]. This technique, probably originating in Persia in the
1400s, presses hydrodynamic correlations into the service of art. Fortunately
for printers, and distressingly for statisticians, a single realization does notresemble the blurry diﬀusion equation.

4.1

Eddy diﬀusion of a front

Figure 5 shows a single realization of the evolution of a “front” under the
RW advection process. The front is the sharp border which separates whitefrom dark; initially this line coincides with the y-axis. We suppose that theconcentration is c = −1 for x < 0 and c = +1 for x > 0. Successive pulsesof the renovating wave produce an increasingly folded front and the c = −1
ﬂuid invades the region x > 0 in long thin tendrils. The central question is:
How well is the process in ﬁgure 5 described by the diﬀusion equation?
We know that given many realizations of this process, the long-time ensemble average of these realizations will follow the diﬀusion equation ct =Dcxx , with the initial conditions c(x, 0) = ±1. The solution of this problem
isc = erfη ,

where

x
η= √ .
2 Dt(25)

Figure 6 shows this smooth erf solution which, of course, looks nothing likeﬁgure 5. If the dark ﬂuid in ﬁgure 5 contained radioactive contaminant, andwe wanted to estimate the maximum exposure of at some value of x > 0,then the erf solution in (25) is not useful.
On the other hand, diﬀusivities are useful if we want to know how many
particles are at such-and-such a distance from their initial location. Thus,
ﬁgure 7 shows a histogram of the positions of 10,000 particles which all starton the line x = 0 (the initial front). The Gaussian curve in ﬁgure 7 is

11

t=3τ

t=6τ

t=9τ

t=12τ

t=15τ

t=18τ

t=21τ

t=24τ

Figure 5: Evolution of a front under the advection by the RW model. The front
initially coincides with the y-axis.

the corresponding prediction for the PDF of positions which is obtained bysolving (19) with the initial condition c = δ(x):


1
x2
τ∗
c(x, t) = √
exp −
,
D= .
(26)
4Dt
8
4πDt
The histograms converge slowly to this Gaussian prediction. This asymptotic
success shows that the diﬀusion equations correctly predicts the dispersionof particles when tτ∗ .
An amusing aspect of the simple problem in ﬁgure 5 is that we can easily
calculate the RMS ﬂuctuations of c around the ensemble average concentra12

10.5

erf(η) and 1-erf 2(η)

0
-0.5
-1
-2

-1.5

-1

-0.5

0

η

0.5

1

1.5

2

Figure 6: Evolution of the ensenble-averaged concentration c and its varianceduring the evolution of the front underthe RRW model. Note how most of thevariance is localised around x = 0.
tion in (25). Because c = ±1 we have c2  = 1. Therefore, deﬁning the
ﬂuctuation asc = c − c ,(27)

c  = c2  − c2 = 1 − erf 2 (η) .

(28)

we have
2

The variance c 2  is also indicated in ﬁgure 6.

4.2

Coarse grained averages and spatial ﬁlters

The process in ﬁgure 5 is translationally invariant in the y-direction and sousing only a single realization we can calculate a spatially averaged concentration L
1c̄(x, t) ≡ limc(x, y, t)dy .
(29)
L→∞ 2L −L
The evolution of c̄ will be asymptotically described by the diﬀusion equation.
In a general case, in which there is no statistical symmetry along a particular direction, one can take a single realization and deﬁne a coarse-grained
or low-pass ﬁltered concentration by:ĉ(x, t) ≡ K(x − x )c(x , t) d2 x ;K(|x|)is a ﬁlter.
(30)
13

t=1

t=2

t=3

t=4
150

150
400

100
100

100
200

0

50

50

-0.5

0

t=5

0.5

0
-2

0

0

2

t=10

50

-2

0

0

2

t=20

200

100

0

2

0

5

0

50

t=25

150
100

150

100

100

50

50

50

50
0

-2

-4

-2

0

t=50

2

0

-5

0

t=100

150
150
100

100

0

10

0

-10

0

-5

0

t=400

0

5

150

150

100

100

50

50

50

50
0
-10

0

5

0

10

-20

0

20

40

0

-5

-50

t=1600Figure 7: At t = 0 the front in ﬁgure 5 is tagged by placing 10000 particles along
the y-axis. The historgram above shows the subsequent x-locations of these marker
particles as the front is distorted by the RW model with τ∗ = 1. The Gaussian
curve is given by (26).

The hope is that scale separation between the width of the erf and the swirls
will ensure that ĉ ≈ c. Thus the kernel of the ﬁlter, K in (30), might be a
Gaussian with a width which is at once much smaller than the thickness ofthe erf transition zone and much greater than an individual swirl in ﬁgure 5.
Scale separation is essential here because the ﬁltering operation deﬁned
by the convolution in (30) is not strictly an “average”. Some of the properties
we take for granted when we use averages arec  = 0 ,c = c ,ab = ab .

For the ensemble average, as indicated in (31), everything works.
14

(31)

For a ﬁlter, such asˆin (30), we can deﬁne the ﬂuctuation concentrationc in analogy with (27):c ≡ c − ĉ .

(32)

But then c = 0 and none of the other desiderata in (31) follow. In otherwords, spatial ﬁltering instead of the ensemble averaging introduces a hostof extra assumptions which should be carefully assessed (but almost never
are).

4.3

A digression: Brownian bugs

I have hinted darkly at problems associated with spatial ﬁlters. These issues
are largely ignored by optimistic scientists. The hope is that scale separation
justiﬁes the application of diﬀusive closures to the coarse-grained version
of a single realization. Perhaps a justiﬁcation of this optimistic approachis that the alternative seems so repellant. Nonetheless, it is important torealize that interpreting coarse-grained distributions as ensemble averages
involves a nontrivial assumption. The best way of exposing this assumption
is to exhibit a problem in which spatial ﬁlters and ensemble averages arevery diﬀerent. Accordingly, as a model of biological processes, we consider
random walkers which both die and reproduce. We refer to these biologicalwalkers as Brownian bugs.
The model is formulated by ﬁrst placing N
1 Brownian bugs randomly
in the unit square; the boundary conditions are periodic in both directions.
Each cycle of the simulation begins with a random walk step in which bug
k, located at xk = (xk , yk ), is displaced to a new position(xk , yk ) = mod [(xk , yk ) + (δxk , δyk ); 1] .

(33)

In (33), δxk and δyk are Gaussian random variables and the “mod” is to enforce the periodic boundary conditions and keep each bug in the unit square.
After this random walk step, the second part of the cycle is a “coin toss”
which results in either death (heads) or division (tails). When a lucky bug
divides, the oﬀspring is placed at the same position as the parent. This cycleof random displacement and random birth/death is repeated many times inorder to simulate many generations of reproduction, death and dispersion.
The simulation shown in ﬁgure 8 was implemented in MATLAB using
these rules. The striking result is that the density of bugs spontaneously
15

develops large-scale clumps and voids. Figure 8 seems to show an inverse
cascade of patch sizes: patches emerge on small scales in panel (b) and then,after more cycles, panels (c) and (d) show that the patches have expandedin scale. To quantify this impression, we have computed one-dimensional
concentration spectra which show that an increasingly red spectrum develops.
A seemingly innocuous ingredient of the brownian-bug model is thatdeaths can occur anywhere, but births are always adjacent to a living bug.
This asymmetry between birth and death is crucial for the spontaneous development of the voids and patches evident in ﬁgure 8: if one simulates birthby randomly placing the new bugs in the unit square then no patches form.
This subtle point shows that making the births coincide with living bugs —surely a realistic feature of the model — introduces pair correlations. From
another perspective, one can view the voids in ﬁgure 8 as the result of random extinctions which create voids. The step length of the random walk inﬁgure 8 is such that diﬀusion is not strong enough to ﬁll in the voids createdby extinction.
The ensemble average of the Brownian bug process is described byct = D∇2 c + (λ − µ)c .

(34)

where λ is the birthrate and µ the deathrate. However if the coin-toss is fairthen births and deaths are equiprobable and consequently λ = µ. In thiscase the solution of (34) which satisﬁes the initial condition is
λ = µ,⇒

c = 1/N .

(35)

The uniform density above is the correct answer for the ensemble average concentration: the location of the voids and patches in ﬁgure 8 are accidently
determined by the MATLAB random number generator. If we ensemble average many such patterns then the patches and voids must disappear because
the process is spatially homogeneous.
On the other hand, the spatial average of a single realization, such as thatin ﬁgure 8, will still show concentration patches2 . Thus, in this Brownian
bug example, ĉ = c. Indeed, the patches are surely an important featureof the “real” answer. The correct but useless result in (35) exposes a failureof ensemble averaging. What do we make of this example? Are biological2

If the width of the kernel, K in (30), is larger than the dimension of the patches thenﬁltering will remove the patches. However, since the patches expand in scale, eventually
they will become so large that they survive the blurring power of the ﬁlter.

16

(a) Initial condition, N=20,000(b) 10 cycles, N=19,692

0.8

0.8

0.6

0.6
Y

1

Y

1

0.4

0.4

0.2

0.2

0

0

0.5
X

0

1

0

(c) 100 cycles, N=17,728

0.5
X

1

(d) 1000 cycles, N=21,814

0.8

0.8

0.6

0.6
Y

1

Y

1

0.4

0.4

0.2

0.2

0

0

0.5
X

0

1

0

0.5
X

1

Figure 8: (a) The initial condition is N = 20, 000 randomly located bugs in theunit square. Panels (b), (c) and (d) then show the development of patches after
10, 100 and 1000 cycles of random displacement followed by random birth/death.
As the panel titles indicate, there are random ﬂuctuations in the total size, N , ofthe population. The RMS step length of the underlying random walk is δx2k 1/2 =δyk2 1/2 = 0.005.
problems, with reproduction and death, so fundamentally diﬀerent from theadvection-diﬀusion of chemical tracers? I am not prepared to answer thatquestion in these lectures and I leave further development of this example tothe students.

17

5

Variance budgets

In this section we return to basics and present an alternative view of eddydiﬀusivity. The following arguements emphasize the importance of the concentration variance equation.

5.1

The Reynolds’ decomposition

We depart from the advection-diﬀusion equation
ct + u·∇c = κ∇2 c + s ,(36)

where κ is the molecular diﬀusivity of c and u is an incompressible (∇·u = 0)velocity ﬁeld. Also in (36) we have included a source term, s(x, t), which
forces the system.
The velocity u in (36) is a single realization selected from an ensemble ofvelocity ﬁelds. Then we can introduce the “Reynolds’ decomposition”:
c = c + c ,(37)

where  is the ensemble average and c is the ﬂuctuation from c whicharises in a single realization. Taking the ensemble average of (36) gives
ct + u·∇c + ∇·u c  = κ∇2 c + s .

(38)(The source s is taken to be deterministic, s = s.)
Subtracting the ensemble average in (38) from (36) gives the ﬂuctuation
equationct + u·∇c + ∇·[u c − u c ] − κ∇2 c = −u ·∇c .

(39)

Equation (39) has been organized by taking the source term to the righthand side. Thus we see that advective distortion of the mean gradient, ∇c,
generates the ﬂuctuation c .

5.2

Consequences of linearity

If c = 0 at t = 0 then, because (39) is linear, c and ∇c will be linearly
related. It follows that the eddy ﬂux u c  will also be linearly related tothe mean gradient ∇c. These simple observations, in alliance with the18scale separation assumption, can be used to extract a surprising amount of
information [2].
The scale separation it is plausible that this linear relation between eddy
ﬂux and mean gradient can be developed in a series of the formu c i = −Dij ∗ c,j − Dijk ∗ ∇c,jk + · · ·
(1)

(2)

(40)

The ∗ in (40) indicates that the product also involves convolutions in time,such as t(1)(1)
Dij ∗c,j = Dij (t )c,ij (t − t ) dt .
(41)
0

If the mean ﬁeld is varying slowly over an eddy decorrelation time then the
convolution above approximates to
 ∞
(1)
(1)
 
u c  ≈ −Dij ∗c,j ≈ −Dij (t ) dt c,j (t) .
(42)0

In the simplest cases3
 ∞(1)
Dij (t ) dt = De δij ,

(43)

0

where De is the eddy diﬀusivity. Using (43) the ﬂux gradient relation isu c  − κ∇c = −D∇c ,D ≡ De + κ ,(44)

and the evolution of the average concentration is determined byct ≈ D∇2 c + s .

5.3

(45)

The Gx-trick

The tensors D(n) (t) are determined by the linear operator on the left handside of (39). Thus, these tensors depend on (i) the statistical properties ofu ; (ii) the mean advection u; (iii) the molecular diﬀusion κ. The essential
3“Simple” means that the velocity ensemble is isotropic, homogeneous and reﬂexionallyinvariant. The last requirement means that the mirror image of a particular realization(1)of u is just as probable as u . If the ensemble is reﬂexionally invariant then Dij is asymmetric tensor. This subtle point will be illustrated later in this lecture series.

19

point is that these tensors do not depend on c. At least for the ﬁrst term
(1)in the series, Dij , we can exemplify this by noting that there is a special
solution of (36) in which u = s = 0 and concentration has the formc = Gx + c .

(46)

In (46) the mean concentration is simply c = Gx and the ﬂuctuation c isdetermined from a reduced version of (39):ct + u ·∇c − κ∇2 c = −Gu .

(47)

As emphasized above, the advection of the mean gradient appears as a source
term for c on the right hand side of (47). Because (47) is linear, and G isconstant, the solution c will be proportional to the large-scale gradient G
and otherwise independent of G.
This Gx-trick enforces the platonic ideal of scale separation between theeddies and the mean ﬁeld. If the concept of an eddy diﬀusivity is to haveany validity, then it must work in the simpliﬁed context of (47). In fact, theGx-trick is used in doubly-periodic turbulence simulation to calculate eddy
diﬀusivities. In that context, u = (u, v) and c are eﬃciently represented by
Fourier series. Then (47) is solved using a spectral code and the eddy ﬂux isestimated by computing the integral

 
−1
u c  = A(48)u c dx dy ,over the computational domain. (In (48) A is the total area of the domainso 1 = 1). Notice that in (48) the ensemble average is identiﬁed withan integral over the domain. Later in these lectures we will use this sameprocedure to analytically calculate the eddy diﬀusivities of some spatiallyperiodic velocity ﬁelds.

5.4

The concentration variance equation

An equation for the concentration variance,
1 2
Z ≡ c  ,(49)
2
is obtained by multiplying (39) by c and ensemble averaging. The result is
1
2
Zt + u·∇Z + ∇ ·  u c  − κ∇2 Z = −κ∇c ·∇c  − u c ·∇c . (50)220The terms on the left hand side of (50) can be interpreted as ﬂuxes of Z. The
two terms on the right hand side of (57) are respectively a source of variancedue to advective distortion of the mean gradient, and dissipation of varianceby molecular diﬀusion κ.

5.5

Heuristic closure arguments

In (50) there are three terms which we would like to relate to the meanquantities c and Z. First, there is −u c  · ∇c = De ∇c · ∇c. The
remaining two terms are u c 2 /2 and κ∇c ·∇c .
The correlation u c 2 /2 in (50) is an eddy-ﬂux of c 2 , just as uc  is an
eddy ﬂux of c . Thus, by analogy with (44), we can argue that1  2
u c  = −De ∇Z .
2

(51)

This heuristic argument is discussed further in appendix B.
The ﬁnal term in (50) is the dissipation of variance by molecular diﬀusivity, κ∇c ·∇c . The simplest closure assumption we can make about thisterm is thatκ∇c ·∇c  ≈ βZ ,(52)

where β has the dimensions of time. The closure above relies on dimensional
analysis and the linearity of (36). However, in anticipation of a later discussion of the Batchelor spectrum, we now make some heuristic arguments insupport of (52) which suggest that β is independent of the molecular diﬀusivity.
Suppose that the mean ﬁeld c has a length scale L and that the velocity
ﬁeld u has a length scale Lu (in the RW example Lu = k −1 ). The scale
separation assumption is that
L

Lu .

(53)

The inequality in (53) is exempliﬁed in idealized case of (46) in which L is
inﬁnite. If follows that advective distortion of ∇c generates c ﬁrst on thescale Lu. Then, following our arguments in lecture 1, the scale of c will be
exponentially reduced, like exp(−γt), where γ is roughly proportional to the21

RMS strain of u . This exponential contraction continues till the cascade ishalted by molecular diﬀusion at the scale

κ
-≡
.
(54)
γ
Using arguments from lecture 1, we can estimate that the time taken for thisarrest at - ist ≈ γ −1 ln (Lu/-) .

(55)

Then the smallest length scale in thec -ﬁeld is - and, plausibly, the gradient√is ∇c ∼ cRM S /- where cRM S ≡ 2Z. We now have a simple estimateκ∇c ·∇c  ∼ γZ. This rough argument leads to the closure in (52), withβ ∝ γ, and the caveat that t > t .
We can summarize the arguments above by rewriting the variance equation (50) as
Zt + u·∇Z − D∇2 Z = De ∇c·∇c − βZ ,(if t ≥ t ) .

(56)

The most dubious approximation is probably (52). To conclude this discussion we will interpret the variance equation in two speciﬁc examples.

5.6

Example 1: the dispersing front

First consider the dispersing front in ﬁgure 5. In this example s = κ = u =0 and we have already know from (28) thatZ=

1
1 − erf 2 (η) ,2x
η= √ .
2 Dt(57)

On the other hand, since κ = 0, it follows that D = De and β = 0. With
these simpliﬁcations the variance equation (50) reduces toZt − DZxx = D∇c·∇c ,(58)

where c is the erf-solution in (25). As a consistency check, one can showthat (57) is the solution of the variance equation in (58).
This example shows that the destruction of variance by molecular diﬀusivity is not required in order to prevent an accumulation of variance: thesource on the right hand side of (58) is balanced by eddy diﬀusion.
22

6

30

5

25

4

20

y

y

35

3

15
2
10
1

5
0

0

10

20
x

0

30

0

2

4

6

x

Figure 9: A simulation of the source problem, with s = cos(x/6), using the RW
model to generate u in (36). There is no molecular diﬀusivity (κ = 0). The
left hand panel shows the whole domain (the length of the side is 12π) while
the right hand panel shows a smaller subdomain (the length of the side is 2π).
The concentration ﬁelds were generated by 10 pulses of the renovating wave usingτ∗ = 3 (that is, t = 30).

5.7

Example 2: a large-scale source

In this second example the tracer is injected by a source s = cos qx in (36).
We also take u = 0 so that the mean concentration ﬁeld is obtained bysolving1 −Dq 2 tct − D∇2 c = cos qx , ⇒ c =1−ecos qx .
(59)
Dq 2(To apply the diﬀusion equation the scale of the source, q −1 , must be muchlarger than the scale of the velocity ﬁeld.) A steady mean concentration
pattern is established when Dq 2 t
1.
The concentration variance is determined by solving the variance equation(56)
2
1 De 
2−Dq 2 tZt − D∇ Z =(1 − cos 2qx) − κ∇c ·∇c  .
(60)
1−e
2
2
2D q
In (60), the solution in (59) has been used to evaluate the source term on
the right hand side and we have left the diﬀusive sink in its exact form.
23

Concentration at t= 200τ=400
300
concentration

200
100
0
-100
-200
-300
0

10

20

30

40

50

60

50

60

x

cRMS at t= 200τ=400
200

RMS

180

c

160
140
120
0

10

20

30

40
x

Figure 10: This ﬁgure compares analytic results with a numerical solution of(62), taking q = 1/10, and using the RW model to generate u. The persistence
parameter is τ∗ = 2 and the results are at t = 400 (that is, 200 renovation cycles).
The concentration c(x, 400) is calculated on a 400 × 400 grid using the method
in appendix B. In the top panel there are three curves: the concentration as afunction of 0 < x < 20π along the line y = 0 (the jagged dotted curve); they-averaged concentration deﬁned in (64); the analytic√ result in (59) (the smoothsinusoid). The bottom panel compares the cRM S = 2Z obtained by solving (60)analytically with cRM S estimated using (64).

It is clear from (60) that the molecular diﬀusion, κ, plays an importantrole. If κ = 0 then the long time solution of (60) has a component which
eventually grows linearly with time:κ = 0,⇒

Z ∝ t/2Dq 2 .
24

(61)

Thus, without molecular diﬀusion, there is “runaway variance”. Ultimately,
in a single realization, the mean ﬁeld in (59) will be buried under enormous
ﬂuctuations.
To give an intuitive derivation of (61) we argue that with κ = 0 theconcentration on each ﬂuid element is determined by solving the Lagrangian
equationDc= cos qx(t) ,
Dt

(62)

where x(t) is the randomly changing x-position of the particle. Thus, theconcentration on each particle is undergoing a random walk along the c-axis,
which is induced by the random motion of the particle through the cos qx
source function. The decorrelation time of this walk is the time it takes aparticle to diﬀuse through a distance of order q −1 , which is 1/Dq 2 . Thus, ina time t, there are roughly N (t) ∼ Dq 2 t independent steps along the c-axis.
But because the source acts coherently for a time 1/Dq 2 with a strength oforder unity, the step length of this random walk is roughly ∆c ∼ 1/Dq 2 .
Thus, the mean square displacement of c is:c  ∼ (∆c)2 N (t) ∼
2

t
,
Dq 2

(63)

which is the ﬁnal result in (61).
It is interesting to compare the analytic results in (59) and (60) with anumerical solution of (62). Thus we must compute the spatial averages


1 L
1 L
2
c̄(x, t) ≡c(x, y, t) dy,cRM S (x, t) ≡[c(x, y, t) − c̄]2 dy , (64)L 0
L 0using the numerical solution, and compare these with the analytic results forc and Z = c2RM S /2. The best way to make this comparison is to obtainc(x, y, t) on a regular grid in the (x, y)-plane. As a bonus, one can then alsouse contouring routines to make pretty pictures of the concentration ﬁeld
(see ﬁgure 9).
The concentration ﬁeld is calculated on a regular grid using the procedure described in Appendix C (essentially the method of characteristics).
Figure 10 shows good agreement between this simulation and analytic results. Notice that in ﬁgure 10 the variance Z peaks where ∇c is greatest.
This illustrates that concentration ﬂuctuations are produced by advective
25

distortion of the mean gradient: where the mean gradient is large there islots of variance. But Z = 0 even where ∇c = 0 (for example, at x = 0 and
x = 10π in ﬁgure 10). Thus, where the source term on the right hand sideof (60) vanishes, the diﬀusive term D∇2 Z supplies variance.

5.8

Cautionary remarks

In the both examples above there is no molecular diﬀusion (κ = 0) andconsequently there is no destruction of variance by the term κ∇c · ∇c  in(50). As a project for a student, include molecular diﬀusion in the RW model(perhaps by pulsing diﬀusion in alternation with advection) and assess theeﬃcacy of this process. In particular, can the closure in (52) be justiﬁed?

6

Calculation of the RW Green’s function

In this appendix we present an alternative calculation of the RW ensemble
averaged Green’s function, g(r), in (17). The unaveraged Green’s function,
G(x, x0 , t), is the solution of (10). Because the process is spatially homogeneous it is harmless to take x0 = 0 so that
G(x, 0, τ∗ ) = δ [x − τ∗ s sin ϕ] δ [y + τ∗ c sin ϕ] ,(65)

where (s, c) ≡ (sin θ, cos θ). The ensemble average of (65) is computed byintegration over ϕ and θ, as in (4). It is very pleasant that there are two
integrals and two δ-functions. Thus, we ﬁrst do the ϕ-integral by noting that
δ [x − τ∗ s sin ϕ] is nonzero at the two values of ϕ where sin ϕ = x/τ∗ s, andat those positions:d(66)[x − τ∗ s sin ϕ] = ± τ∗2 s2 − x2 .
dϕUsing the standard rule for changing variables in a δ-function, we ﬁnd thatthe average of (65) over ϕ alone isGϕ =1 δ(y + cot θx).
π τ∗2 sin2 θ − x2

(67)

The second integral over θ is performed by noting that δ(y+cot θx) is nonzero
at the two values of θ where cot θ = −y/x, and at those positions
sin2 θ =x2,
x2 + y 2dx2 + y 2[y + x cot θ] = −
.
dθ
x
26

(68)

After changing variables in the δ-function we recover g(r) in (17).

7

Eddy diﬀusion of variance

Ignoring small molecular diﬀusion (κ = 0), if c satisﬁes the advection equation then any function of c satisﬁes the same equation. That is to say
Dc
= 0,
Dt

⇒

Df
= 0,
Dt

(69)

where f (c) = c2 , or exp(c) etcetera. Taking an ensemble average, and making
the same arguments for f (c) as for c, we have thatf t = D∇2 f  .

(70)

In the particular case f = c2 /2, f  = c2 /2 + Z and (70) reduces to
Zt = D∇2 Z + D∇c·∇c .

(71)

Matching the terms in (71) with those in (50) we conclude that u c 2 /2 =
−D∇Z.

8

Numerical simulation of the RW process

Drawing ﬁgures 9 and 10 requires that we obtain the solution of (62) on aregular grid in the (x, y)-plane. This is an opportunity to use the method ofcharacteristics and learn some MATLAB programming techniques.
Equation (3) shows how the movement of a particle in the RW velocity ﬁeld is equivalent to a random map which determines the position at(n + 1)τ∗ in terms of the previous position at nτ∗ . If this particle carries
a concentration, c(x, t), which changes because of the cos qx source in (62),then the concentration changes can also be calculated and expressed as amap in discrete time.
Thus, suppose that the concentration on a particle at time t = nτ∗ is cn .
Then the change in concentration during nτ∗ < t < (n + 1)τ∗ is obtained byintegratingDc= cos [qxn + qun (t − nτ∗ )] ,
Dt
27

(72)

where the x-velocity of the particle is un = sn sin(cn xn + sn yn + ϕn ), with(sn , cn ) ≡ (sin θn , cos θn ) is the constant velocity of the particle. The integralof (72) can be written as
cn+1 = cn +sin(qxn+1 ) − sin(qxn )
.
qun

(73)

With equations (73) and (3) we can integrate forward in time and so determine the concentration on a particle at t = nτ∗ .
However we need to determine the concetration at t = nτ∗ at a speciﬁed
grid point x. So the trick is to...
%% Solution of%%%%
Dc/Dt=cos(qx);%%%% cos(q x) is a large -scale source and u is the RW velocity.
%% The RW streamfunction is psi=cos[cos(theta) x+ sin(theta)y + phi]clc
N=400;%% Use an N*N grid in the plotting window
q=1/6;%% The wavenumber of the cos q x source
LL=2*pi/q;%% LL is the domain size
npulse=10%% The number of renovation cycles
tau=3;%% The pulse duration of the wave%% Lwin is the side of the square plotting window.
%% Set Lwin=LL to see the big picture. To see small scale details,%% try Lwin = 2*pi. We draw two subplots with different Lwin’s
nloop=0;
for Lwin=[LL 2*pi]
nloop=nloop+1
x=linspace(0,Lwin,N);%% x is the coordiante in the plotting window.
h=x(2);%% The grid spacing in the plotting window
for j=1:N
jj=[((j-1)*N+1):(j*N)];pos(jj,1)=zeros(N,1)+(j-1)*h;pos(jj,2)=x’;
28

end
conc=zeros(N*N,1);%% The position of the N^2 particles are stored in pos with%% N^2 rows and 2 columns. Each vertical segment of%% length N in pos contains particles with the same initial x-position.
%% the column vector conc contains the concentration on the%% N*N particles in pos. Initially, conc=0 at the N*N%% grid points. Then we integrate%% backwards in time to find the concentration change.
for k=1:1:npulse
theta=rand*2*pi;
wavevec=[cos(theta),sin(theta)]’;phase=rand*2*pi;
vel=sin(pos*wavevec+phase)*[wavevec(2),-wavevec(1)];conc=conc-sin(q*pos(:,1))./(q*vel(:,1));
pos=pos+tau*vel;
conc=conc +sin(q*pos(:,1))./(q*vel(:,1));end%% Emerging from this loop, we have the the new positions%% and the new concentrationconc=reshape(conc,N,N); %% conc is reshaped into an N*N matrix
hh=subplot(1,2,nloop)
colormap(’gray’)imagesc(x,x,conc)
axis equal
xlabel(’x’)
ylabel(’y’)
axis([0 Lwin 0 Lwin])set(hh,’ydir’,’norm’)end29References[1] S. Childress and A. D. Gilbert. Stretch, Twist, Fold: The Fast Dynamo.
Springer, Berlin, 1995.
[2] H.K.Moﬀatt. Transport eﬀects associated with turbulence with particularattention to the inﬂuence of helicity. Rep. Prog. Phys., 46:621–664, 1983.
[3] W.A. Medeiros. Marbling Techniques. Watson-Guptill, New York, 1994.
[4] G.I. Taylor. Eddy motion in the atmosphere. Phil. Trans. Roy. Soc.
London, 215:1–16, 1915.

30