1

Chapter 2

Diffusion Equation
Part 1
Environmental Transport and Fate
Benoit Cushman-Roisin
Thayer School of Engineering
Dartmouth College
Recall: Flux of a substance
The flux of a substance in a particular direction
is defined as the quantity of that substance
passing through a section perpendicular to that
direction per unit area and per unit time:
t A
m
q
A
=

= =
duration time area sectional - cross
section - cross through passes t amount tha
flux
If the amount is a mass, then the flux q is a rate of mass per area per time,
and its units are kg/m
2
s.
2
For the moment (first part of Chapter 2), we restrict our attention to the diffusive
component of the flux. Thus,
j u c q
q
+ =
+ = diffusion advection
In general, the flux consists of two components:
advection (= passive movement by carrying fluid)
diffusion (= random movement with respect to mean motion of fluid)
due either to molecular agitation
or turbulent fluctuations in fluid flow.
j q =
only.
We saw that the diffusive component of the flux can be expressed as
dx
dc
D j =
in which D is the diffusivity coefficient
and dc/dx is the concentration gradient along the length of the domain.
We now use this form of flux into the mass budget for a short portion [x, x+Ax] of a
one-dimensional domain.
A t x x q A t x q
dt
dc
V
) , ( ) , (
Export Import
A + =
=
Mass budget
Since the volume V of the slice is equal to AAx, we can write after division by V:
x
t x q t x x q
dt
dc
A
A +
=
) , ( ) , (
And, in the limit of a vanishingly short distance Ax:
x
q
t
c
c
c
=
c
c
[A switch from total (straight) to partial (curved) derivative is warranted
because there are now more than one variable in the derivatives at this stage.]
3
Putting the pieces together
With the earlier formulation of the diffusive flux
x
c
D q
c
c
=
the budget equation becomes
x
q
t
c
c
c
=
c
c
|
.
|

\
|
c
c
c
c
+ =
c
c
x
c
D
x t
c
This equation is the 1D diffusion equation. It is occasionally called Ficks second law.
In many problems, we may consider the diffusivity coefficient D as a constant.
In that case, the equation can be simplified to
2
2
x
c
D
t
c
c
c
=
c
c
Remember, however, that in an environmental context D is likely to represent
the effect of turbulent fluctuations, which may have spatially varying statistics.
(Ficks first law)
The preceding equation can hardly be solved in any general way. Each solution depends
critically on boundary and initial conditions, specific to the problem at hand.
2
2
x
c
D
t
c
c
c
=
c
c
First and foremost, we need to know how many initial and boundary conditions are
necessary so that the problem is neither underspecified or overspecified.
The time derivative is of first order. This demands that we provide one and only one
condition in time. Typically, the initial concentration distribution is given:
0 at ) (
0
= = t x c c
The spatial derivative on the right-hand side is of second order and thus calls for two
boundary conditions. Typically, any problem will include one boundary condition at each
end of the domain (say x = x
1
and x = x
2
). While it would be easier from a mathematical
perspective to have the end concentrations imposed, such as
2 2
1 1
at ) (
at ) (
x x t c c
x x t c c
= =
= =
it is much more typical to encounter flux boundary conditions.
4
Flux boundary conditions are:
2 2
1 1
at ) (
at ) (
x x t q
x
c
D
x x t q
x
c
D
= =
c
c

= =
c
c

A common situation is that of an impermeable boundary.

In that case, there cannot be any flux across the boundary,
and there is no concentration gradient at that boundary:
0 =
c
c
x
c
at an impermeable boundary
Mixed conditions (e.g., concentration given at one end of the domain and
flux specified at the other) are possible.
Prototypical 1D solution
The diffusion equation is a linear one, and a solution can, therefore, be
obtained by adding several other solutions. An elementary solution (building block)
that is particularly useful is the solution to an instantaneous, localized release
in an infinite domain initially free of the substance.
Mathematically, the problem is stated as follows:
Infinite domain: < x < +
D = constant
No initial concentration, except for a localized release: c(x) = M o(x) at t = 0
Since the substance will take an infinite time to reach the infinitely far ends of the
domain, we impose:
0 lim lim = =
+
c c
x x
at finite times (t < ).
In the above, Mis the total mass of the substance released per unit cross-sectional area,
and o(x) is the Dirac function [o(x) = 0 for x 0, o(x) = + at x = 0, with area under the
infinitely tall and infinitely narrow peak being unity].
5
We expect the substance to spread gradually from its location of release and
the concentration distribution therefore to look like this:
After some rather lengthy mathematics (see written notes), we find that the solution is:
|
|
.
|

\
|
=
Dt
x
Dt
M
t x c
4
exp
4
) , (
2
t
We note that this solution already meets the boundary conditions (vanishing
concentrations far away on both sides). The initial condition is also met away from x = 0,
because when t = 0, the exponent is and the exponential vanishes.
Conservation of the total amount of the substance is what sets the front coefficient.
It requires that
M dx t x c dx t x c = =
} }
+

+

) , ( ) , (
0
And one can show that the above solution satisfies this requirement.
6
Random-walk analogy
Random-walk process
Take a one-dimensional domain (1D axis), divide it in a series of boxes of identical
lengths Ax (bins), and place one particle in one of the bins. Imagine now that this
particle is endowed with a mechanism that makes it jump randomly every time interval,
At, according to the following rules:
There is a 25% chance that the particle will hop one bin to the left,
There is a 25% chance that it will hop one bin to the right, and
There is a 50% chance that it will remain in the same bin.
This is the 1D version of the more general 2D random walk, commonly referred to as the
drunkard's path: Late into the night, an inebriated fellow leaves a bar and has no
recollection of where he is and where he is going; every second or so, he makes a step
forward, backward (he turns completely around), leftward or rightward, making a random
path that may look like that displayed in the figure below.
The question is:
Where is this drunkard expected to be after m steps?
7
Because of the random nature of the problem, the answer can only be in terms
of probabilities. We thus define
p(nAx, mAt)
as the probability that the particle is in bin number n at time t = m At.
We can calculate the probability at one time step based on the previous time step:
At time (m + 1) At, the probability that the particle is in bin n is equal to the probability
that it was already there [ p(nAx, mAt) ] times the probability that it stayed there (50%),
plus the probability that it was one bin to the left [ p((n-1)Ax, mAt) ] times the probability
that it jumped to the right (25%), plus the probability that it was one bin to the right
[ p((n+1)Ax, mAt) ] times the probability that it jumped to the left (25%).
Mathematically, we have:
] , ) 1 [(
4
1

] , ) 1 [(
4
1

) , (
2
1
] ) 1 ( , [
t m x n p
t m x n p
t m x n p t m x n p
A A + +
A A +
A A = A + A
And, the initial condition is:
0 for 0 ) 0 , (
1 ) 0 , 0 (
= = A
=
n x n p
p
The solution to this discretized system is:
0.0029 0.016 0.054 0.121 0.193 0.226 0.193 0.121 0.054 0.016 0.0029 m=6
0.0010 0.010 0.044 0.117 0.205 0.246 0.205 0.117 0.044 0.010 0.0010 m=5
0.004 0.031 0.109 0.219 0.273 0.219 0.109 0.031 0.004 m=4
0.016 0.094 0.234 0.313 0.234 0.094 0.016 m=3
0.063 0.250 0.375 0.250 0.063 m=2
0.250 0.500 0.250 m=1
1.000 m=0
+5 +4 +3 +2 +1 0 -1 -2 -3 -4 -5 n=
Note the spreading over time and the thinning of the fringes.
8
To compare with the continuous solution obtained from the diffusion equation, we sample
the continuous solution at x = nAx and t = mAt :
|
|
.
|

\
|
A
A

A
= A A
t Dm
x n
t Dm
M
t m x n c
4
) (
exp
4
) , (
2
t
Then, let us set (Ax)
2
equal to 4DAt to clean up the exponent:
|
|
.
|

\
|

A
= A A
m
n
x m
M
t m x n c
2
exp ) , (
t
Finally, let us set the front coefficient to match the probability value at
(n = 0 , m=5):
|
|
.
|

\
|
= A A
m
n
m
t m x n c
2
exp
5
246 . 0 ) , (
The c values thus calculated compare quite favorably to the probabilities of the
random-walk process, at least for the larger times.
0.0029 0.016 0.054 0.121 0.193 0.226 0.193 0.121 0.054 0.016 0.0029 m=6
0.0010 0.010 0.044 0.117 0.205 0.246 0.205 0.117 0.044 0.010 0.0010 m=5
0.004 0.031 0.109 0.219 0.273 0.219 0.109 0.031 0.004 m=4
0.016 0.094 0.234 0.313 0.234 0.094 0.016 m=3
0.063 0.250 0.375 0.250 0.063 m=2
0.250 0.500 0.250 m=1
1.000 m=0
+5 +4 +3 +2 +1 0 -1 -2 -3 -4 -5 n=
0.0035 0.016 0.050 0.115 0.190 0.225 0.190 0.115 0.050 0.016 0.0035 m=6
0.0017 0.010 0.041 0.111 0.201 0.246 0.201 0.111 0.041 0.010 0.0017 m=5
0.005 0.029 0.101 0.214 0.275 0.214 0.101 0.029 0.005 m=4
0.016 0.084 0.228 0.318 0.228 0.084 0.016 m=3
0.053 0.236 0.389 0.236 0.053 m=2
0.202 0.550 0.202 m=1
m=0
+5 +4 +3 +2 +1 0 -1 -2 -3 -4 -5 n=
p values
c values
We conclude that the random-walk and diffusion are similar processes (Einstein, 1905).
= anchor between
both representations
9
c(x)
x
A graphical iteration
(An old pencil-and-paper method predating computers or even slide rules!)
A curve c(x) is given graphically as a
succession of points, Ax apart, in the (x, c)
plane (heavy dots), and line segments
connecting adjacent points approximate
the curve (solid line).
The following graphical constructions are then made:
1. Join the midpoints of adjacent line segments (red lines),
2. At each location, move the function value to the level of the red line (purple lines),
3. Repeat these steps again and again.
The obvious result of this manipulation is the cutting of spikes, filling of valleys and overall
smoothing of the curve.
Ax
At every cycle of this procedure, we effectively replace the current value of the function by
a new value, which is the average of the averages of the neighboring values. Thus,
) (
4
1
) (
2
1
) (
4
1

2
) ( ) (
2
) ( ) (
2
1
) ( new
x x f x f x x f
x x f x f x f x x f
x f
A + + + A =
(

A + +
+
+ A
=
In other words, we are performing a (, , ) running average, which is the same
as following the random walk.
We conclude:
Diffusion, random-walk and curve smoothing are all similar.
Diffusion random walk diffusion spreads.
Diffusion curve smoothing diffusion smoothes.
10
Spreading
Spreading induced by diffusion implies that the spatial extent of contamination grows with
time. For practical purposes, we wish therefore to quantify this spreading, viz. to have an
answer to the question:
How wide is the zone affected by the contaminant?
We easily conceive that the quantity that tells us how wide is a diffusing patch of
contamination should be a function of time. However, what should be its precise
definition is not trivial. Indeed, the prototypical solution
|
|
.
|

\
|
=
Dt
x
Dt
M
t x c
4
exp
4
) , (
2
t
exhibits spreading over time but does not yield a specific width because c is nonzero all
the way to infinity, starting immediately after the moment of release. [This instantaneous
infinitely wide effect is obviously unrealistic and is attributed to the mathematical simplification of a
continuum medium.] Therefore, looking for the edges of the concentration distribution does
not yield an answer. We need to do something else.
To determine objectively the width of a pollutant patch, we then resort to integral
quantities. Let us explore a few of them.
A first integral is the area under the curve:
constant
-
= =
}
+

M dx c
giving the total amount of substance; this tells nothing about the width of the patch.
A second integral is formed by inserting the coordinate x, to inject the notion of distance,
0
-
=
}
+

dx xc
which vanishes by symmetry; this tells us where the mean position of the patch is
( for this solution) but still nothing about the width.
A third integral is
0 = x
0 ) (
-
2
=
}
+

dx c x x
11
0 ) (
-
2
>
}
+

dx c x x
is a function of time, positive and always nonzero.
Because represents the squared distance to the mean position, the above
integral says something about the average distance to the center of the patch and
therefore holds information about the patch width.
The time evolution of this width should then provide information about the rate of
spreading.
0 ) (
1
-
2 2
= =
}
+

dx c x x
M
o
2
) ( x x
dx c M
}
+

=
-
in which
and dx xc
M
x
}
+

=
-
1
For convenience, we define the normalized quantity
For the prototypical solution above, calculations yield
Dt
Dt
x
d e
Dt
dx
Dt
x
Dt
M
x
M
2
4
with
4

4
exp
4
1
2
2
2
2 2
=
= =
|
|
.
|

\
|
=
}
}
+

+

, , ,
t
t
o
,
Thus, Dt 2 = o
2
t
=
We immediately note that the distance o
grows in time, but as the square root of
time rather than time itself. Thus, the
spreading goes on with time, never
stopping but gradually slowing down.
12
The next question is by which numerical factor (such as 4, 1/2, or whatever)
should this o be multiplied to yield a practical definition of the patch width.
The c(x,t) distribution at any given time obeys the so-called Gaussian (or bell) curve.
Normalizing distance x by o and concentration c by its maximum, central value
the function becomes
Dt
M
c
t 4
max
=
(
(

|
.
|

\
|
=
2
max
2
1
exp
o
x
c
c
which contains no parameter. Graphically, the curve is universal.
It is the so-called Gaussian bell curve encountered in probability and statistics.
2.5%
2.5%
95%
o o 2 2 + < < x
The region contains 95% of the pollutant, which is about most of it.
We can then say that the overall width of the region over which the pollutant is spread
is
Dt
Dt
66 . 5
2 4 4 width
=
= = o
We call this the 4-sigma rule.
13
A remark
There are occasions when one should distinguish between
4o the full width, and
2o the spreading distance from the location of release.
Take for example a smokestack plume.
A pertinent question is:
By which distance does the ground level
become affected?
The answer is: where the underside of
the plume begins to touch the ground,
that is, where 2o = H, the height of the
stack.
ut x
H Dt H
=
= =

2 2 2o
D
uH
x
8
2
=
An example
A tank aboard a barge traveling along the Chicago Ship Canal suddenly collapses,
releasing its benzene content (C
6
H
6
, density = 0.879 g/cm
3
), of which 100 liters find their
way quickly to the water. The rest of the benzene remains contained on the barge.
Assuming rapid mixing across the canal section (8.07 m deep and 48.8 m wide) and
estimating the turbulent diffusion coefficient at 3.0 m
2
/s, what are the concentrations of
benzene 2, 6, 12 and 24 hours after the accident, at the site of the spill and 300 m away?
To solve this problem, we first determine the mass of benzene that was spilled. Since
the density of benzene is 0.879 g/cm
3
= 0.879 kg/L, this mass m is:
Over the cross-section of the canal, we have
kg 9 . 87 L 100 kg/L 879 . 0 volume density mass = = = = m
2
m
kg
0.2232
m) m)(48.8 (8.07
kg 9 . 87
area sectional cros
mass
= = =

=
A
m
s
M
14
The concentration at the location of the spill (x = 0) is given by
) hrs in (
mg/L 606 . 0
) sec in (
kg/m 0364 . 0

4
) ( ) , 0 ( ) (
3
max
t t
Dt
M
t c t c t c
spill
= =
= = =
t
0.124 mg/L 24 hours
0.175 mg/L 12 hours
0.247 mg/L 6 hours
0.428 mg/L 2 hours
Concentration
at spill
Time
since spill
The peak concentration decreases over time. This can be expected since diffusion
spreads the pollutant away from the point of release.
At a distance x = 300 m away, the concentration is obtained from the full formula
|
|
.
|

\
|
=
Dt
x
Dt
M
t x c
4
exp
4
) , (
2
t
with
) hrs in (
083 . 2
) s / m 0 . 3 ( 4
m) 300 (
4
2
2 2
t t Dt
x
= =
Thus
|
|
.
|

\
|
=
) hrs in (
083 . 2
exp
) hrs in (
mg/L 606 . 0
) , m 300 (
t t
t c
0.113 mg/L 24 hours
0.147 mg/L 12 hours
0.175 mg/L 6 hours
0.151 mg/L 2 hours
Concentration
300 m away
Time
since spill
Note how the concentration at 300 m first
increases (because it takes time for the
benzene to diffuse over that distance) and
then decreases (because further diffusion
leads to dilution).
15
Observe how, at
some distance
from the peak, the
concentration first
increases and then
decreases.
We can actually calculate the exact time at which the concentration reaches its
maximum and what that maximum is.
Setting to zero the time derivative of solution, we obtain the time t
max
of the maximum
concentration for any distance x from the release location:
D
x
t
2
2
max
=
Substitution in the solution then yields the
maximum concentration:
x
M
x
M
D
x
t x c x c
2420 . 0
2
1
exp
2 2
, ) (
2
max
=
|
.
|

\
|
=
|
|
.
|

\
|
= =
t
In the case of the benzene spill in the Chicago Ship Canal,
with M= 0.2232 kg/m
2
and x = 300 m:
c
max
= 0.180 mg/L at t
max
= 4 hrs 10 min.