kt = 25 m2

= 0.8

= 2.5

0.003

0.0005

(c) (f)

2

2

2

2

0.0

0.0

5 0 5 10 15

0 5 10 15

x (m)

x (m)

In many applications, diffusion is one com-

Fig 2.1 Evolution of a topographic scarp, illustrating (a)

ponent of a more complex model. The transport

elevation, (b) slope, and (c) curvature. In (a), arrows of

of contaminants in groundwater ¬‚ow, for ex-

varying length represent the sediment ¬‚ux at each point. In

the diffusion model, the ¬‚ux is proportional to the local ample, can be quanti¬ed using a combination

slope, and the resulting raising or lowering rate of the surface of advection (i.e. the translocation of contami-

is proportional to the change in ¬‚ux per unit length, which, in nants as they are carried along with the ¬‚uid)

turn, is proportional to the curvature. (d) and (e) Graphs of

and diffusion (i.e. the spreading of contaminant

elevation, slope, and curvature for ¬ve times following scarp

concentration from hydrodynamic dispersion as

offset (κt = 0.01, 0.03, 0.1, 0.3, 1.0.)

the contaminant travels along variable ¬‚ow paths

through a porous aquifer). Similarly, the concen-

tration of dust downwind from a playa is the re-

weathering-limited slopes and channels evolve

sult of both advection and turbulent diffusion.

advectively. Analytic and numerical techniques

Advection describes the motion of the dust as it

for solving the advection equation will be pre-

is carried downwind, while diffusion describes

sented in Chapter 4.

2.1 INTRODUCTION 33

the spreading of the plume as it is mixed by at-

mospheric turbulence. In many of the examples

in this chapter, we will assume a ¬xed chan-

nel head at the base of a diffusing hillslope. Of

course, channels have their own complex evolu-

nonlinear

q

tion and are rarely ¬xed at the same elevation transport

over time. As such, a complete understanding of

hillslopes or channels requires a fully coupled

model of the ¬‚uvial system of which hillslope

diffusion is one component. Although the tech- linear (slope-proportional)

transport (i.e. diffusion eqn.)

niques used in this chapter are framed around

solving a single equation that has limited ap-

Sc

plication in geomorphology, the techniques we

will discuss form the basis for solving many real-

world problems in which diffusion is one part of Fig 2.2 Schematic graph of ¬‚ux versus slope. The diffusion

a more realistic model. equation assumes a linear relationship between ¬‚ux and

The applicability of the diffusion model to slope. More generally, sediment ¬‚ux increases nonlinearly

hillslope evolution depends on the processes act- with slope, diverging at a critical slope value Sc . However,

even the most general ¬‚ux relationship can be approximated

ing to move sediment on the hillslope. Creep

as linear for small slopes. As such, the predictions of linear

and rain splash are examples of hillslope pro-

and nonlinear transport models are the same for

cesses that are accurately represented by the dif-

gently-sloping topography.

fusion equation for gentle slope angles (Carson

and Kirkby, 1972). In these processes, particle

movement takes place by a series of small jumps

the divide (Carson and Kirkby, 1972):

triggered by freeze/thaw cycles, wetting/drying,

or momentum imparted from rain drops. In these ‚h

q = ’κ x (2.6)

jumps, particles may move both upslope and ‚x

downslope but the distance traveled by a parti-

cle moving downslope is slightly greater than the The evolution of hillslopes governed by this

distance of those traveling upslope. Moreover, the ¬‚ux relationship is considered in Section 2.3.2.

increase in distance traveled by a particle mov- Mass movements are not diffusive because the

ing downslope increases linearly with the hill- sediment ¬‚ux increases rapidly and nonlinearly

slope angle or gradient for relatively small slopes. with small increases in hillslope gradient near

the point of slope failure Sc , as illustrated in

Bioturbation is another important hillslope pro-

Figure 2.2.

cess that can often be approximated as diffusive,

The value of κ controls how fast diffusion

particularly for gravelly hillslopes of arid envi-

takes place. In the Basin and Range Province of

ronments where processes such as slope wash

the western US, an approximate value of κ =

are relatively ineffective due to the coarse tex-

1 m2 /kyr has often been used based on studies of

ture of the slope material. Slope wash, rilling,

pluvial shorelines and other landforms of known

and mass movements are examples of hillslope

age (Hanks, 2000). It is widely recognized, how-

processes that are not diffusive. Slope wash is

ever, that the 1 m2 /kyr value varies spatially ac-

not diffusive because the shear stress exerted by

cording to climate, vegetation, soil texture, and

overland ¬‚ow increases with distance from the di-

other variables, even if that variation is not yet

vide as ¬‚ow accumulates along the slope. While

well characterized. One drawback of the diffusion

hillslopes dominated by slope wash are not gov-

approach in geomorphology is that we often do

erned by the diffusion equation, they can be de-

not know the value of κ for a given location with

scribed by a more general, spatially variable dif-

any precision. This is a valid criticism, but it is

fusion equation in which the ¬‚ux is proportional

still important to study hillslope evolution even if

to the hillslope gradient and the distance from

34 THE DIFFUSION EQUATION

the absolute rates of change cannot be precisely uplifting mountain range. In such cases, the

determined. diffusion equation simpli¬es to the time-

Three classic publications form the basis of independent equation:

many of the results in this chapter. Carslaw and ‚ 2h

κ +U = 0 (2.7)

Jaeger (1959) is the de¬nitive text on analytic so- ‚ x2

lutions to the diffusion equation. Although the

where U is the uplift rate. Solving for the hill-

book was written primarily for engineering ap-

slope curvature, Eq. 2.7 gives:

plications to heat conduction in solids, the solu-

‚ 2h U

tions presented in that book can be applied to =’ (2.8)

‚ x2 κ

diffusion problems that arise in any ¬eld of sci-

ence. Culling (1963) was among the ¬rst authors Integrating Eq. (2.8) gives an equation for the hill-

to apply the diffusion equation to hillslope evo- slope gradient:

lution. His paper is still the most complete ref- ‚h U

= ’ x + c1 (2.9)

erence on the subject. Culling™s paper includes ‚x κ

solutions appropriate for many different initial

where c 1 is an integration constant. In order to

conditions as well as a stochastic treatment of

constrain the value of c 1 , it is necessary to spec-

the relationship between the random motion of

ify the value of ‚h/‚ x at one end of the hillslope

sediment particles on the hillslope and the diffu-

pro¬le. The upslope end of the hillslope is cho-

sive evolution of the hillslope. Carson and Kirkby

sen to coincide with a divide in this case. Divides

(1972) is the most complete reference on the rela-

are defined as locations across which no water or

tionship between speci¬c hillslope processes and

sediment passes. As such, a no-flux boundary con-

their signature forms.

dition is appropriate for the upslope end of this

In the following sections we present solutions

pro¬le. The no-¬‚ux boundary condition is given

to the diffusion equation primarily within the

by

context of hillslope evolution. As such, the dif-

‚h

fusing variable h will represent elevation along =0 (2.10)

‚x

a topographic pro¬le. It should be emphasized, x=0

which implies c 1 = 0. Integrating Eq. (2.9) gives

however, that the solutions we describe are gen-

eral. Therefore, any solution that we obtain for h

U2

h(x) = ’ x + c2 (2.11)

(for a particular initial condition and boundary

2κ

conditions) can be applied to any other diffusion

where c 2 is another integration constant. At the

problem, whether the physical process is heat

downslope end of the pro¬le, the elevation is as-

conduction, contaminant transport, etc. That is

sumed to be constant. This boundary condition is

part of the power of applied mathematics: once

appropriate for a hillslope-channel boundary in

a set of solutions or a solution method has been

which all of the sediment delivered to the chan-

learned it can often be rapidly applied to many

nel is transported away from the slope base with

different physical systems.

no erosion or deposition. If the hillslope has a

length L and an elevation of zero at the slope

base:

2.2 Analytic methods and

h(L ) = 0 (2.12)

applications

then c 2 = U L /2κ and the hillslope pro¬le is given

by

2.2.1 Steady-state hillslopes

U2

Steady-state landscapes have received a great deal h(x) = (L ’ x 2 ) (2.13)

2κ

of attention in geomorphology in recent years.

The topography of a steady-state landform is time- Equation (2.13) shows that diffusive hillslope pro-

independent, such as when uplift and erosion ¬les are parabolas in cases of steady-state up-

are in perfect balance at every point in a rapidly lift. The hillslope relief in Eq. (2.13) is given by

2.2 ANALYTIC METHODS AND APPLICATIONS 35

U L 2/2κ, so the relief of diffusive hillslopes is pro- Fourier showed that nearly any function could be

represented as a series of sinusoidal functions:

portional to uplift rate and slope length, and

inversely proportional to diffusivity. Equations ∞

nπ x

f (x) =

(2.9)--(2.13) illustrate how boundary conditions an sin (2.16)

L

are used to solve differential equations. n=1

Equation (2.8) suggests that in regions of sim- where f (x) is any function that goes to zero

ilar climate and hillslope processes, hillslope cur- at x = 0 and x = L . Given a function f (x) that

vature can be a relative measure of tectonic up- obeys these boundary conditions, the values of

lift rates. Equation (2.8) may be applicable to the the Fourier coef¬cients an can be calculated using

early stage of uplift when hillslope gradients are

L

2 nπ x

still relatively low and weathering rates are suf- an = f (x ) sin dx (2.17)

L L

¬ciently rapid for hillslopes to be regolith cov- 0

ered. In general, however, Eq. (2.8) is of limited The solution to the diffusion equation with ini-

use in many areas of rapid uplift because the re- tial condition given by Eq. (2.16) is the summa-

quirements for the diffusion equation generally tion of Eq. (2.15) from n = 1 to n = ∞:

do not apply. Later in this chapter we will solve

∞

nπ x ’κn2 π 2 t/L 2

the solution to the nonlinear diffusion equation.

h(x, t) = an sin e (2.18)

L

That equation usually provides a more accurate n=1

representation of hillslope evolution in rapidly

Substituting Eq. (2.17) into Eq. (2.18) gives the gen-

uplifting terrain.

eral solution

∞

L

2 nπ x

2.2.2 Fourier series method h(x, t) = f (x ) sin

L L

Two general classes of solutions exist for the dif- 0 n=1

nπ x ’ κn2 π 2 t

fusion equation: series solutions and similarity

— sin e L2 dx (2.19)

solutions. First consider the one-dimensional (1D) L

diffusion equation in a region bounded by x = 0

This approach suggests a powerful means of

and x = L . By 1D, we mean a pro¬le that varies

solving the diffusion equation in bounded re-

with only one independent spatial variable, x.

gions: using the initial condition f (x), solve for

Let™s assume that the initial condition at t = 0

the coef¬cients an of the Fourier series repre-

is given by a function f (x) and that the value of

sentation of f (x); then, plug the an values into

h is set to zero at both ends of the region for

Eq. (2.18) to obtain the solution. Other boundary

all time t (i.e. h(0, t) = h(L , t) = 0). If the initial

conditions (e.g. values of h(0, t) and h(L , t) that

condition f (x) is given by

are nonzero, or boundary conditions speci¬ed by

nπ x hillslope gradients) can also be handled by using

f (x) = an sin (2.14)

L cosine, linear, and/or constant terms in the sum-

mation. For example, if we wish to place a divide

then the solution to the diffusion equation is

at x = 0 then the series solution is

nπ x ’κn2 π 2 t/L 2

h(x, t) = an sin ∞

e (2.15) L

1 nπ x

h(x, t) =

L f (x ) cos

L 2L

’L n=1

The correctness of this solution can be checked

nπ x ’ κn2 π22 t

— cos

by differentiation and substitution of the re- e 4L dx (2.20)

2L

sults into the diffusion equation. This solution

is only valid for the sinusoidal function given by Equation (2.20) is the same as Eq. (2.19) ex-

Eq. (2.14), so it may not appear to be of much use. cept that cosine terms are used and the hill-

slope domain is assumed to be from x = ’L

Equation (2.15), however, is actually a powerful

to x = L . If f (x) is chosen to be symmetric (i.e.

building block that can be used to solve the dif-

f (x) = f (’x)) then Eq. (2.20) is consistent with

fusion equation for nearly any initial condition

by using a Fourier series representation for f (x). the divide boundary condition given by Eq. (2.13).

36 THE DIFFUSION EQUATION

The Fourier series approach is only applicable integration. The approach is to introduce a new

to bounded regions. Of course, all geomorphic variable that is a combination of x and t:

applications are bounded, because the Earth is x

·= √ (2.22)

2 κt

not in¬nite in extent! In some cases, however, it

is useful to de¬ne mathematical models in in¬- and rewrite the diffusion equation for θ and its

nite or semi-in¬nite domains. For example, if one boundary conditions in terms of ·. This can be

is interested in modeling the ¬rst few years of done (following the approach in Turcotte and

radionuclide migration into an alluvial deposit Schubert (1992)) using the chain rule for differ-

that is hundreds of meters thick, a semi-in¬nite entiation:

model is more convenient to work with and prac-

‚θ ‚θ ‚· ‚θ 1·

1x1 dθ

= = ’√ = ’

tically speaking more accurate than a Fourier

‚t ‚· ‚t ‚· 4 κt t d· 2t

series approach.

(2.23)

The Fourier series approach works because

‚θ dθ ‚· dθ 1

= = √

the diffusion equation is linear (i.e. it does not (2.24)

‚x d· ‚ x d· 2 κt

involve squared or higher-order terms of h or

‚ 2θ 1 d2 θ ‚· 1 1 d2 θ

its derivatives). Only in linear equations can dif- =√ = (2.25)

‚ x2 2 κt d·2 ‚ x 4 κt d·2

ferent solutions be superposed to obtain other

The diffusion equation for θ becomes

solutions.

1 d2 θ

dθ

’· = (2.26)

2.2.3 Similarity method 2 d·2

d·

Series solutions with ¬xed-elevation boundary

and the boundary conditions become θ(∞) = 0

conditions are appropriate for hillslopes with a

and θ (0) = 1. Equation (2.26) is a second-order or-

well-de¬ned base-level control (such as a nearby

dinary differential equation (ODE) that can be

channel that carries sediment away from the

reduced to a ¬rst-order ODE by introducing the

slope base to maintain a constant elevation). In

variable

some cases, however (such as when sediment is

dθ

allowed to deposit and prograde at the slope base) φ= (2.27)

d·

it is more appropriate to assume that diffusion

occurs over an in¬nite or semi-in¬nite region. In Equation (2.26) then becomes

such cases, similarity solutions are a powerful 1 dφ

’·d· = (2.28)

approach. Let™s consider the semi-in¬nite region 2φ

given by x = 0 to x = ∞ with an initial condition

Equation (2.28) can be integrated to obtain

h(x, 0) = h 0 . Let™s also assume that at t = 0 the el-

’·2 = ln φ ’ ln c 1

evation h is instantaneously lowered to zero at (2.29)

x = 0. In the context of hillslope evolution, this

Taking the exponential of both sides of Eq. (2.29)

problem corresponds to instantaneous base level

gives

drop of an initially planar hillslope or terrace.

dθ

In order to solve this problem, it is useful to 2

φ = c 1 e’· = (2.30)

introduce a dimensionless elevation given by d·

Integrating Eq. (2.30) gives

h

θ= ’1 (2.21) ·

h0 2

e’· d· + 1

θ = c1 (2.31)

The diffusion equation for θ is identical to the 0

where · is an integration variable and the con-

diffusion equation for h. The boundary condi-

tions on θ are now given by θ(x, 0) = 0, θ(0, t) = 1, dition θ (0) = 1 was used to constrain the integra-

and θ(∞, t) = 0. Similarity solutions make use of tion constant. Application of θ(∞) = 0 requires

a mathematical trick that reduces the diffusion that

∞

equation (a partial differential equation) to an 2

e’· d· + 1

0 = c1 (2.32)

ordinary differential equation that is solved by 0

2.2 ANALYTIC METHODS AND APPLICATIONS 37

√

The integral in Eq. (2.36) is given by π/2, so the on h. The initial condition for w is given by

√

constant c 1 has the value ’2/ π. Equation (2.31)

U2

w(x, 0) = ’ (L ’ x 2 )

becomes (2.39)

2κ

∞

2 2