. 4
( 14)


numerical model results
u = 5 m/s, K = 5 m2/s, p = 0.05 m/s
wind direction
(d) (e)
numerical model results
numerical model results
u = 5 m/s, K = 5 m2/s, p = 0.05 m/s
1 km
u = 5 m/s, K = 5 m2/s, p = 0.05 m/s

q = 10
q = ’10 q = 30
10 20 40
0 80 cm

Fig 2.19 (a) Plot of eolian silt thickness versus downwind distance, with analytic solutions for the two-dimensional model for
representative values of the model parameters. (b) Schematic diagram of model geometry. Depositional topography shown in this
example is an inclined plane located downwind of source (but model can accept any downwind topography). (c) Color maps of
three-dimensional model results, illustrating the role of variable downwind topography. In each case, model parameters are
u = 5 m/s, K = 5 m2 /s, and p = 0.05 m/s. Downwind topography is, from left to right, a ¬‚at plane, an inclined plane, and a
triangular ridge. Width and depth of model domain are both 6 km. (d) Color maps of three-dimensional model results for
deposition downwind of Franklin Lake Playa, illustrating the role of variable wind direction. From left to right, wind direction is
θ = ’10—¦ (0—¦ is due south), 10—¦ , and 30—¦ . (e) Map of three-dimensional model results obtained by integrating model results over a
range of wind conditions, weighted by the wind-rose data in Figure 2.18a. For color version, see plate section. Modi¬ed from
Pelletier and Cook (2005).

topography because mass loss from the plume is the spatially variable deposition that occurs as
still assumed to occur along a horizontal plane in the plume intersects with complex topography.
Eq. (2.68). However, by using the elevated plume Two-dimensional (2D) model solutions iden-
concentration -- i.e., c[x, y, h(x, y)] -- to estimate tify a fundamental length scale for downwind de-
the deposition term pc, the model approximates position. Smith (2003) showed that Eq. (2.68) can

be approximated by a power-law function along model corresponding to the wind-rose diagram
the plume centerline: of Figure 2.18a. The thickness values mapped in
Figure 2.19e correspond to a uniform Q value of
Q x
c(x, 0, 0) = √ 1.0 m/kyr (assuming a Qa3 surface age of 50 ka),
2 π pL p Lp
obtained by scaling the model results for differ-
ent values of Q to match the maximum thickness
where L p is a characteristic length scale given by
L p = uK / p2 . Figure 2.19a, for example, presents values between the model and observed data.
The modeled and observed deposition patterns
two plots of Eq. (2.70) with different sets of model
match each other in most respects. However, the
parameters. These parameter sets were chosen
secondary hot spot located on the western pied-
to have different values of K and p but simi-
mont of the study area is not reproduced by
lar values of L p . The similarity of the resulting
the model. This feature could result from more
plots indicates that the pattern of 2D dust depo-
westerly paleowinds transporting more dust east-
sition is not sensitive to particular parameters;
ward from the Amargosa River relative to modern
it is only sensitive to the scaling parameter L p .
conditions. Alternatively, it could result from en-
Therefore, raising the value of K and p simulta-
hanced hillslope erosion of Qa3 surfaces in this
neously, as in Figure 2.19a, leads to solutions that
are equally consistent with observations down-
wind of Franklin Lake Playa. This nonuniqueness
makes model calibration (in 2D or 3D) dif¬cult.
2.3 Numerical techniques and
The model is best calibrated using wind-speed
data to constrain u and K to the greatest ex-
tent possible, because these values can be in-
dependently determined using readily-available
2.3.1 Forward-Time-Centered-Space
data. The deposition velocity, in contrast, can-
not be readily determined. By constraining u and
Analytic solutions are elegant but they are often
K , however, model results can then be run for a
only loosely applicable to real problems. Linear
range of values of p to determine the value most
fault scarps and volcanic cones are special cases
consistent with observations.
where the landform symmetry allows for an ana-
Model results are shown in Figure 2.19c to il-
lytic approach because a 2D landform is reduced
lustrate the role of complex downwind topogra-
to 1D. Numerical solutions are generally needed
phy in controlling deposition patterns. In each
in order to model the evolution of landforms that
case, a circular source with uniform Q was con-
sidered with model parameters u = 5 m/s, K = lack symmetry or are controlled by complex base-
5 m2 /s, and p = 0.05 m/s. From left to right, the level boundary conditions.
In order to solve partial differential equations
downwind topography was considered to be ¬‚at,
on a computer, we must discretize the equa-
an inclined plane, and a triangular ridge. Rela-
tions in both space and time. Discretization in
tive to ¬‚at topography, the inclined plane leads
space means that we represent the spatial do-
to a short plume, and the triangular ridge diverts
main as a discrete set of grid points (e.g. from
and bisects the plume.
position x1 to xn for 1D problems). Discretization
Long-term dust transport does not correspond
in time means that the solution will be solved for-
to a single wind direction as assumed in the sim-
ward in discrete intervals of time. First, we must
pli¬ed cases of Figure 2.19c. The effects of multi-
determine how to best represent ¬rst, second,
ple wind directions can be incorporated into the
and higher derivatives in space within this dis-
model by integrating a series of model runs with
crete framework. The starting point is the Taylor
a range of wind directions, each weighted by the
frequency of that wind direction in the wind-
rose diagram. Figure 2.19d illustrates model re-
‚h ( x)2 ‚ 2 h
sults for the Eagle Mountain region with wind h(x + x, t) = h(x, t) + +
‚x 2 ‚ x2
directions from θ = ’10—¦ to +30—¦ (θ = 0—¦ is due x,t x,t

+ O (( x) ) 2
South). Figure 2.19e illustrates the integrated (2.71)

The simplest formula for the ¬rst derivative is grid-point values at the previous time step. The
obtained by neglecting second-order (( x)2 ) terms FTCS method is limited by the fact that it is sta-
and rearranging Eq. (2.71) to obtain the forward- ble only for small timesteps.
difference approximation: The FTCS method starts with an initial value
h i everywhere on the grid. In addition to the ini-
h(x + x, t) ’ h(x, t) ‚h tial condition, boundary conditions on the value
= + O ( x) (2.72)
x of h or its ¬rst derivative must also be speci¬ed

at both ends of the grid. These boundary con-
This is called forward difference because we
ditions may be constant or may vary as a func-
look forward to x + x to compute the differ-
tion of time. Equation (2.76) is then applied se-
ence. Equation (2.72) is not unique, however, be-
quentially from one end of the grid to the other
cause we can also calculate the ¬rst derivative
to specify the grid of values at the next time
using a backward difference:
h(x, t) ’ h(x ’ ‚h The FTCS method is generally not useful for
x, t)
= + O ( x) (2.73)
‚x large grids or higher dimensions because very
x x,t
small step sizes must be taken in order to main-
We can also use a centered-difference approxima- tain stability. To see this, consider the sequence of
tion: grid points in Figure 2.20. Initially, the values of
h are zero everywhere on the grid. At time zero,
h(x + x, t) ’ h(x ’ ‚h
x, t)
= + O (( x)2 ) the left-side boundary condition on h is raised to
2x x,t
a ¬nite value. The right-side boundary condition
is h n = 0. During the ¬rst time step, all of the
grid points remain zero except for the second
The second derivative can be calculated by using
grid point from the left. The second derivative
a forward difference to calculate the slope in the
of h, (h 0 ’ 2h 0 + h 0 )/( x)2 is positive, which re-
positive direction and subtracting the backward 0 1 2
sults in a ¬nite value for h 1 . Note that if the time
difference in the negative direction: 1
step is large enough, the value of h 1 can actually
h(x + x, t) ’ h(x, t) h(x, t) ’ h(x ’
1 x, t) 0
rise higher than that of h 0 . If this happens, the

x x x model will become unstable and will continue
h(x + x, t) ’ 2h(x, t) + h(x ’ x, t) to diverge from the correct solution in subse-
( x)2 quent timesteps. In order to maintain stability,
the timestep t must be less than ( x)2 /2D . The
‚ 2h
= + O (( x)2 ) (2.75)
‚ x x,t
2 time required for diffusion to propagate a dis-
tance L is ≈ L 2 /D . Therefore, the number of time
The discretization of the diffusion equation, steps required to model diffusion is ≈ L 2 /( x)2 .
using Eq. (2.75) and introducing abbreviated dis- To model the diffusion equation in 1D on a grid
crete notation h(i x, n t) = h in , is given by with 1000 grid points, therefore, requires a min-
imum of about one million time steps. Once we
h i+1 ’ 2h in + h i’1
h in+1 ’ h in n n
=D move to 2D and 3D problems with tens or hun-
( x)2
dreds of thousands of grid points, one million
= h in + (h ’ 2h in + h i’1 )
h in+1 n
(2.76) time steps becomes computationally unfeasible
( x)2 i+1
on even the fastest computers. Another related
This Forward-Time-Centered-Space (FTCS) method disadvantage of the FTCS method is that infor-
is a useful means for solving the 1D diffusion mation travels very slowly from one end of the
equation for small grids. This scheme is numer- grid to the other. Figure 2.20 illustrates that it
ically stable provided that the time step is less would take 1000 time steps for the grid points
than or equal to ( x)2 /2D . Equation (2.76) is also on the right side of the grid to feel any effect
known as an explicit scheme because the value of the boundary condition on the left side of
of each grid point is an explicit function of the the grid. This slow propagation of information

h h
unstable, ”t too large


implicit method:
explicit method:

Fig 2.20 Illustration of the difference between (a) explicit
Eq. (2.45)). Substituting Eq. (2.45) into Eq. (2.2)
numerical schemes and (b) implicit schemes. In explicit
schemes, grid values at the previous time steps are used to
solve for values at the next time step. In the example, the
left-side boundary, h 0 , is abruptly raised at t = 0. In the ‚h ‚ 2h ‚h
=κ x 2 + (2.77)
explicit scheme, only the ¬rst grid point to the right, h 1 will
‚t ‚x ‚x
feel the effects of this change at the boundary in the following
time step. If the value of the time step is chosen to be too
Here we consider the slope-wash-dominated ver-
large, the explicit scheme can overcorrect and make h 1 larger
than h 1 . In the implicit scheme, the effect of raising the left sion of instantaneous base-level drop of an ini-
boundary is simultaneously felt throughout the grid no matter tially ¬‚at terrace considered in Section 2.2.5.
how large it is, resulting in a more stable solution. Figure 2.21a compares hillslope pro¬les ob-
tained for the nonuniform diffusion equation
(Eq. (2.77)) with those of the classic diffusion
is one reason why the time step must remain equation (Eq. (2.20)). For early times, the two
so low. models show no difference. For later times, how-
So, if the FTCS is so poor, why use it at all? ever, hillslope erosion occurs more rapidly in the
The answer is that it is easy to program and, in nonuniform diffusion case, but otherwise there
some ways, more versatile than other schemes we is little difference in the pro¬le shape.
will discuss. In the next two sections we utilize These results are broadly consistent with
the FTCS method to solve for the evolution of those of Carson and Kirkby (1972), who compared
hillslopes with spatially variable diffusivity and hillslope pro¬les for uniform and nonuniform
nonlinear sediment transport. cases. Figure 2.21b is based on the classic ¬g-
ure from their book, illustrating the signature
hillslope forms for different hillslope processes.
2.3.2 Evolution of hillslopes with In this ¬gure, Carson and Kirkby argued that
spatially-varying diffusivity slopes governed by the uniform diffusion equa-
tion (m = 0 and n = 1 in their terminology) were
On hillslopes where slope wash is the predomi-
nant transport process, the classic diffusion equa- characteristically convex downward. In contrast,
the slope-wash-dominated case (i.e. m = 1 and
tion cannot be used. In such cases it is appropri-
n = 0) was characterized by a straight pro¬le with
ate to model hillslope sediment ¬‚ux as increas-
ing linearly with distance from the divide (i.e. no curvature.

kt =
20 m2

kt = 0.6
0.6 h
1280 m2 h0
h0 kt =
1280 m2

kt =
kt =
q= k 0.2
0.2 2560 m2
2560 m2 kx
q= L
0 50
0 50
x (m)
x (m)
(a) (b)
slope wash
q= L
without gulying

q= k
soil creep

slope wash
with gullying

0 1.0
where n is a constant. Substituting Eq. (2.78) into
Fig 2.21 (a) Comparison of hillslope evolution of an initially
Eq. (2.2) gives
horizontal terrace with instantaneous base-level drop for
⎡ ¤
homogeneous diffusion (i.e. ¬‚ux is given by Eq. (2.2)) and
‚h ‚h/‚ x n|‚h/‚ x|n
⎢ ⎥
inhomogeneous diffusion (i.e. ¬‚ux is also proportional to
=κ ⎣1 + ¦
n n
distance from the divide, or Eq. (2.45)). The two evolution
1 ’ |‚h/‚ x| ‚h/‚ x|
Scn 1 ’
Sc Sc
equations yield almost identical hillslope forms, but the rate
of erosion is higher in the inhomogeneous case if the same κ (2.79)
value is used in both cases. (b) Characteristic slope pro¬les
Andrews and Hanks (1985) proposed using Eq.
based on the characteristic forms of Carson and Kirkby
(2.79) with n = 1, while Andrews and Bucknam
(1972). These authors argued that uniform and nonuniform
(1987), Roering et al. (1999), and Mattson and
diffusion result in distinctly different hillslope forms, in
Bruhn (2001) assumed n = 2. Appendix 1 provides
contrast to the model results in (a).
a program to solve for hillslope evolution gov-
erned by Eq. (2.79).
2.3.3 Evolution of hillslopes
Figure 2.22 presents the nonlinear hillslope
with landsliding
response (Eq. (2.79)) to instantaneous base-level
The literature on scarp analysis has included a
drop with an initially planar terrace. The ini-
vigorous debate on linear versus nonlinear diffu-
tial phase of the model is an instantaneous re-
sion models. For example, Andrews and Bucknam
laxation of the hillslope to the critical angle
(1987) ¬rst proposed that the ¬‚ux of material is
Sc at time t = 0. During the ¬rst 1 kyr of the
proportional to the gradient for small values, but
model (i.e. κt = 1 m2 assuming κ = 1 m2 /s), little
increases towards in¬nity as the gradient reaches
hillslope rounding takes place and the hillslope
a critical value Sc :
wears back to a slightly lower angle. As time pro-
ceeds, downwearing gives way to diffusive round-
‚h/‚ x
q = ’κ ing of the hillslope crest. The rate of hillslope
|‚h/‚ x|
1’ evolution decreases rapidly with the overall relief

evolution is very different from that of the linear
kt = 0
kt = 0 diffusion equation. Appendix 1 provides a pro-
1 m2 gram to solve for hillslope evolution governed by
0.8 10 m2
Eq. (2.79). This program may be used as a tem-
100 m2
plate to solve hillslope evolution problems with a
h broad range of ¬‚ux equations, initial conditions,
h0 and boundary conditions.
kt = 103 m2
2.3.4 2D Evolution of alluvial-fan terraces
Figure 2.23 illustrates solutions to the 2D diffu-
sion equation using the FTCS method. The model
20 50
10 40
0 domain is 100 m on each side and the hillslope
x (m)
is assumed to diffuse with κ = 1 m2 /kyr. Figure
2.23a illustrates the channel mask grid that iden-
Fig 2.22 Plots of nonlinear hillslope pro¬les (Eq. (2.79)) at
times κt = 1 m2 to 103 m2 following instantaneous base-level ti¬es which grid points are to be kept at a ¬xed
drop of an initially planar terrace. elevation. At t = 0 a planar alluvial terrace is
abruptly incised according to the channel mask
grid. This example is designed to be a 2D analog
because of the strongly nonlinear relationship be- of the abrupt terrace entrenchment considered
tween ¬‚ux and hillslope gradient. This slowing in Section 2.2.5. Figures 2.23b--2.23d illustrate the
is evident in the very large range of time scales terrace morphology after (b) 50 kyr, (c) 120 kyr,
(varying over three orders of magnitude) between and (d) 500 kyr following entrenchment. The mo-
the ¬rst and the last pro¬le shown. Clearly this tivation for this example is the observation that

Fig 2.23 Solution to the diffusion
equation with κ = 1 m2 /kyr in the
neighborhood of a series of gullies
(shown in (a)) kept at constant base
level and a model domain of
0.01 km2 . This model represents the
evolution of an alluvial-fan terrace
abruptly entrenched at time t = 0.
After (b) 50 kyr, diffusional rounding
of the terrace near gullies has

penetrated ≈ κt or 7 m into the
terrace tread and planar terrace
t = 50 kyr
rill “mask”
treads are still widely preserved.
After (c) 120 kyr approximately
(c) (d)
11 m of rounding has taken place.
Finally, after (d) 500 kyr, erosional
processes have removed all planar
terrace remnants and a rolling
ridge-and-ravine topography

t = 500 kyr
t = 120 kyr

latest Pleistocene alluvial fan terraces in the general matrix equation, because the matrix to
western US often have extensive preservation of be inverted in Eq. (2.83) is a tridiagonal matrix.
planar terraces. As terraces increase in age to Tridiagonal matrices have zero values everywhere
several hundred thousand years or longer, they except for on the main diagonal and one row
eventually lose all planarity and become ridge- above and below it. Such matrices can be inverted
and-ravine topography. Appendix 1 provides a very quickly using a short-cut version of Gaussian
program for modeling the evolution of terraces elimination (Press et al., 1992). This approach to
using the FTCS method to obtain the results in solving diffusion problems is often called an im-
Figure 2.23. plicit method, because the unknown grid values
at time n + 1 appear on both sides of Eq. (2.80)
2.3.5 Implicit method rather than just on one side as in explicit
Fortunately, there are ways to solve diffusion methods.
The computational advantage of the implicit
problems much more ef¬ciently than with the
FTCS method. For example, we can rewrite the method is hard to overstate. In the explicit FTCS
right side of Eq. (2.76) in terms of the new time scheme, the new values for each grid point are
step n + 1 rather than the previous time step n: calculated based on the old (out-dated) informa-
tion from the previous time step. This is the rea-
hi = hi + ’ 2h i + h i’1
n+1 n+1 n+1 n+1
h (2.80) son why it takes 1000 time steps for the right
( x)2 i+1
side of the grid to know anything about the left
At ¬rst glance Eq. (2.80) may seem nonsensical: side. In the implicit method, however, all of the
grid values at time n + 1 appear on both sides of new values everywhere on the grid are calculated
the equation. So, how can we solve for h in+1 on the simultaneously. Figure 2.20b, for example, illus-
left side of Eq. (2.80) if it has to be input on the trates how an instantaneous jump in the bound-
right side ¬rst? The answer is that Eq. (2.80) can ary condition at one end of the grid is felt im-
be rearranged to form a matrix equation: mediately at the next time step throughout the
⎛ ⎞
grid. In the implicit scheme, grid points instantly
1 0 0 0
⎜ ’± (1 + 2±) ⎟
’± receive information from every other grid point
⎜ ⎟
⎜ ⎟
’± (1 + 2±) ’± ...⎟
⎜0 no matter how large the grid is. This difference
⎜ ⎟
⎝0 ⎠
’± (1 + 2±) makes the implicit method far more stable than
... explicit methods. In fact, the implicit method is
⎛ n+1 ⎞ ⎛ n ⎞ stable for any size time step. The accuracy of the
h0 h0
⎜ h n+1 ⎟ ⎜ h n ⎟ solution, however, still depends on the time step.
⎜ 1 ⎟ ⎜ 1⎟
⎜ n+1 ⎟ ⎜ n ⎟
(2.81) In practice, it is useful to run two versions of the
— ⎜ h2 ⎟ = ⎜ h2 ⎟
⎜ n+1 ⎟ ⎜ n ⎟ same implicit simulation with time steps that dif-
⎝ h3 ⎠ ⎝ h3 ⎠
fer by a factor of 2. The difference between the
... ...
two solutions provides an estimate of the accu-
where ± = κ t/( x) . Equation (2.81) is a matrix racy of the results for the larger time step.

equation of the form The implicit method can be made more accu-
(2.82) rate by averaging the second derivative at time
Ahn+1 = hn
n + 1 and time n:
which has the solution
h in+1 = h in + h i+1 ’ 2h in+1 + h i’1
n+1 n+1
=A h
n+1 n
h (2.83) 2( x)2
+ h i+1 ’ 2h in + h i’1
n n
The order of the matrix in Eq. (2.83) is the
same as the number of grid points in the model.
Solving the diffusion equation in 1D using Eq. This is known as the Crank--Nicholson method
(2.83), therefore, requires inverting a 1000 — 1000 and it is a work horse for solving diffusion prob-
matrix. That may sound ominous, but, in fact, lems. The implicit method was ¬rst used in a
solving Eq. (2.83) is much easier than solving a landform evolution context by Fagherazzi et al.

(2002), but its use is still surprisingly limited As a concrete example of the ADI technique,
given the power of the method. we consider a model for hillslope evolution in
Hanaupah Canyon. In this model, we input the
USGS 30 m/pixel DEM of Hanaupah Canyon as
2.3.6 Alternating-Direction-Implicit
our starting condition. We then solve the 2D dif-
fusion equation with κ = 10 m2 /kyr on all pix-
The implicit method is particularly powerful in
els where the contributing area is less than a
2D and 3D because the advantage of the implicit
threshold value of 0.1 km’2 . The results of this
scheme grows with the number of grid points
model are shown in Figure 2.24 for a simulation
and grid sizes grow geometrically with the num-
of 10 Myr duration using time steps of 1 Myr. The
ber of dimensions. In 2D, the implicit method
ADI technique is particularly useful for solving
solves the equation
2D and 3D problems with complicated bound-
ary conditions, as in this case, where each chan-
h i, j = h i, j + h i+1, j + h i, j+1
n+1 n+1 n+1
( x) nel pixel is a ¬xed-elevation boundary. It should
’4h i, j + h i’1, j + h i, j’1
n+1 n+1 n+1
be noted that this example is a bit contrived be-
cause the hillslopes in Hanaupah Canyon are, as
Unfortunately, Eq. (2.85) does not yield a tridiag- noted previously, dominated by mass movements
onal matrix because grid points are dependent and therefore are not well represented by the
on points to their right and left but also up and diffusion equation. Nevertheless, there are many
down. The up and down linkages in Eq. (2.85) re- other areas of low relief where this approach
sult in a much more complex matrix compared would be a good model for hillslope evolution.
to the simple tridiagonal matrix. Code for implementing this model is presented in
The power of the implicit method can still Appendix 1.
be brought to bear on this problem, however,
by solving the problem repeatedly as a series
of 1D implicit problems. Each row of the grid
is calculated as a 1D implicit problem with the 2.1 Many glaciers around the world have responded
grid values in the up and down directions (h i, j+1 to global warming by increasing their rates of
n+1 basal sliding. Assuming that basal sliding occurs
and h i, j’1 ) treated as constants. Then, the direc-
as a result of warmer atmospheric temperatures
tions are reversed: each column is calculated as
conducted through the ice, estimate the lag time
a 1D implicit problem with the grid values in
between atmospheric warming and glacier re-
the right and left directions treated as constants
sponse for a glacier 100 m thick.
using the output from the implicit calculations 2.2 Consider a long alluvial ridge of height 100 m,
done previously. This process is repeated several width 1 km, and a sinusoidal cross-sectional to-
times, solving for rows and then columns, until pographic pro¬le (i.e. h(x) = h 0 sin(π x/»)). Assum-
a prescribed level of accuracy has been achieved ing a diffusivity of κ =1 m2 /kyr, what will be the
as indicated by very small improvements to the height of the ridge 1 Myr in the future?
calculated values. One additional advantage of 2.3 As time passes, fault and pluvial shoreline scarps
become more subdued features of the landscape.
this Alternating-Direction-Implicit (ADI) scheme
Given diffusivity values prevalent in the western
is that it can handle complex boundary con-
US of κ ≈ 1 m2 /kyr, fault scarps 10 kyr in age are
ditions almost effortlessly. If, for example, we
readily apparent. Older scarps, however, eventu-
wished to solve the diffusion equation in the
ally become too diffused to identify. Assuming an
neighborhood of a meandering river, we could in-
initial angle of 30—¦ , how old is the oldest identi¬-
put a mask of grid values that are kept constant,
able scarp in the western US? Assume that scarps
assuming that the river acts as a constant base can be identi¬ed if their slope values are greater
level for the surrounding hillslope. Then, the ADI than 2% (relative to the far-¬eld slope).
technique can be used and the ¬xed elevations of 2.4 You have been contracted to conduct a prelimi-
the river can simply be ¬‚agged so that the value nary seismic hazard assessment for an area that
at n + 1 is kept equal to the value at n. has fault scarps of unknown age with a range of

Fig 2.24 Example application of
(a) the ADI technique to modeling
hillslope evolution in Hanaupah
Canyon. (a) Shaded relief of US
Geological Survey DEM of
Hanaupah Canyon. (b) Solution with
10 Myr of hillslope evolution
assuming κ = 10 m2 /kyr and
channels with ¬xed elevations.
Channels were de¬ned to be all
pixels with a contributing area
greater than 0.1 km2 .


heights. You wish to determine which scarps are any initial condition. Begin by constructing two
likely Holocene in age. Develop an expression for columns for x and h ( x, 0). Choose at least ten
the minimum midpoint angle required for a scarp grid points and use any values you wish for x
to be Holocene in age, based on the height of the and h(x, 0) (i.e. you may choose to construct a
scarp, 2a, and the regional diffusivity κ. hypothetical fault scarp). Construct a third col-
2.5 The fallout of man-made radionuclides from at- umn in your spreadsheet for h(x, t), where the
mospheric nuclear testing in the 1950s and 1960s value of t is chosen to be small enough to en-
provided a pulse of tracer particles that scientists sure stability. Assume that the value of h is con-
have used to quantify transport rates in many dif- stant at the two ends of the domain. By copying
ferent settings. Although ocean currents are very the h(x, t) column nine times and modifying the
complex, the vertical distribution of tritium in links accordingly, plot the solutions for h(x, 2 t),
the ocean is approximately diffusive. Plot the rel- h(x, 3 t) . . . h(x, 10 t).
ative concentration of tritium in the ocean as a 2.8 Older moraines often have a lag of coarse material
function of depth assuming that a pulse of tri- on their crests. Develop a model for this coarsen-
ing based on the diffusion equation with a κ value
tium was deposited in the ocean in 1960 and
κ = 100 m2 /yr. that depends on mean grain diameter. For exam-
2.6 In this chapter we used Fourier series methods in ple, assume that the moraine parent material is
1D only, but they can also be used in 2D for sim- 30.
ple geometries. Use the Fourier series method to 2.9 Identify a radiometrically dated cinder cone
solve for the shape of a rectangular ridge of initial in the western US using the published litera-
height H , length L , and width W as a function ture. Download USGS Seamless 10 m/pixel DEM
of time. Plot the solutions in contour or shaded (http://seamless.usgs.gov) for an area that includes
relief for several different time values assuming the cone and extract radial topographic pro¬les of
κ = 1 m2 /kyr. the cone. Determine a local κ value for the cone
2.7 Using a spreadsheet program, implement the using the Fourier--Bessel solutions to determine a
best-¬t pro¬le using different values of κ.
FTCS scheme for the diffusion equation for

2.10 Deserti¬cation is lowering water levels in many of base level be felt 10 and 100 yr into the
lakes worldwide (e.g. Dead Sea). As lakes act as future?
the base level for channels draining into them, 2.11 Using available DEM data (e.g. http://seamless.
lake-level drops will likely trigger channel inci- usgs.gov) extract a topographic pro¬le across the
sion on the margins of these lakes. Model the evo- Sierra Nevada Mountains. Compute the tempera-
lution (analytically or numerically) of a channel ture pro¬le (analytically or numerically) at a shal-
longitudinal pro¬le with diffusivity κ =10 m2 /yr low depth h in the crust beneath the pro¬le. As-
sume a geothermal gradient of 20 —¦ C/km and a
subject to continuous base-level drop of 1 m/yr.
mean annual surface temperature of 0 —¦ C.
How far from the shoreline will the effects
Chapter 3

Flow routing

promise between realism, computational speed,
3.1 Introduction and ease of use.
Flow-routing algorithms have application to
a wide range of geomorphic problems. Most
Digital Elevation Models (DEMs) play an impor-
landscape-evolution models, for example, assume
tant role in modern research in Earth surface pro-
that erosion or sediment transport is a function
cesses. First, DEMs provide a baseline dataset for
of the discharge or drainage area routed through
quantifying landscape morphology. Second, they
each pixel. As such, most landscape evolution
enable us to model the pathways of mass and en-
models require a ¬‚ow-routing algorithm in order
ergy transport through the landscape by hillslope
to quantify erosion and sediment transport. Flow-
and ¬‚uvial processes. Given the importance of
routing algorithms are also important for mod-
DEMs in a broad range of geoscienti¬c research,
eling contaminant transport in ¬‚uvial systems.
the ability to digitally process DEMs should be a
In such cases, ¬‚ow-routing algorithms are used
part of every geoscientists™ toolkit.
to transport contaminants from source regions
Flow-routing algorithms lie at the heart of
and mix them with uncontaminated sediments
DEM analysis. Flow-routing algorithms are used
downstream. In this chapter we describe the com-
to model the pathways of mass and energy
mon ¬‚ow-routing algorithms used in DEM anal-
through the landscape. There is no unique ¬‚ow-
ysis and present three detailed case studies of
routing method because different constituents
their application.
move through the landscape in different ways.
Water, for example, moves through the landscape
somewhat differently than sediment. Also, ¬‚ow-
routing models are necessarily simpli¬ed mod-
3.2 Algorithms
els of transport. Indeed, the only ideal model for
water ¬‚ow in the landscape is the full solution
3.2.1 Single-direction algorithms
to the Navier--Stokes equations of ¬‚uid dynamics.
The simplest ¬‚ow-routing algorithms are single-
Modeling the Navier--Stokes equations in a com-
direction algorithms. These algorithms assume
plex topographic environment, however, is be-
that all incoming ¬‚ow is routed in the direc-
yond the scope of any computer. As such, all ¬‚ow-
tion of steepest descent. By ˜˜¬‚ow,™™ we usually
routing methods involve some approximation to
mean discharge or some proxy for discharge such
the real physics of mass and energy transport.
as contributing area. These models are general,
The key, then, is to develop a suite of ¬‚ow-routing
however, and may apply to any constituent that
algorithms that partition mass and energy
moves through the hillslope and ¬‚uvial system.
down-slope in ways that mimic the complex pro-
The D8 algorithm (O™Callahan and Mark,
cesses of actual ¬‚uid ¬‚ow. Modelers can then
1984) works by searching the eight neighbors
choose which algorithm represents the best com-

(including diagonals) for the steepest down-slope because ¬‚uvial topography is characterized by
gradient. All ¬‚ow from the central pixel is then tributary ¬‚ow in many areas at scales larger than
delivered to that pixel. To use this algorithm 90 m. Today, however, high-resolution DEMs con-
to compute the contributing area at each point structed by LIDAR and low-elevation photogram-
within a grid, each pixel is ¬rst initialized to have metry commonly have resolutions of 1 m/pixel.
a contributing area equal to ( x)2 , where x is Accurate ¬‚ow routing in these high-resolution
the pixel width. Then, moving through the grid DEMs requires multiple-direction ¬‚ow-routing al-
in a systematic manner (e.g. upper-left to lower- gorithms that partition ¬‚ow from a central point
right) the path of steepest descent is computed into multiple down-slope pixels.
from each pixel to the DEM boundary, increment-
ing the value of the contributing area by ( x)2 in
3.2.2 Multiple-direction algorithms
each pixel along each path.
The Multiple-Flow-Direction (MFD) algorithm of
One problem with the D8 algorithm is that
Freeman (1991) was the ¬rst algorithm to par-
¬‚ow pathways predicted by this algorithm have
tition ¬‚ow into multiple down-slope pixels. In
segments that are constrained to be multiples
this algorithm, ¬‚ow from a central pixel is par-
of 45—¦ . Since ¬‚ow pathways in nature can take
titioned into all down-slope pixels, weighted by
on any orientation, at a minimum a good ¬‚ow-
slope. To implement this technique, the pixels
routing algorithm should not depend on the rect-
must be processed in a speci¬c order. In this and
angular nature of the underlying DEM. In prac-
most other multiple-¬‚ow-direction algorithms,
tice, the D8 algorithm tends to produce overly
pixels must be processed from the highest eleva-
straight ¬‚ow pathways since the local hillslope
tions to the lowest. Since mass and energy move
or channel aspect (i.e. orientation) must deviate
downhill, processing pixels from highest to low-
by at least 45—¦ in order for the D8 algorithm to
est elevation ensures that all of the ¬‚ow coming
resolve a change in ¬‚ow direction.
into a pixel from upstream has been computed
The ρ8 algorithm (Fair¬eld and Laymarie,
before the ¬‚ow is partitioned down-slope. To im-
1991) attempts to overcome the problem of overly
plement this algorithm, an index list must ¬rst
straight ¬‚ow pathways in the D8 algorithm.
be created that ranks the DEM pixels from high-
The ρ8 algorithm randomly assigns ¬‚ow from
est to lowest. This ranked list can be computed
the center grid cell to one of its down-slope
very ef¬ciently using the quicksort algorithm
neighbors (including diagonals) with a probabil-
(Press et al., 1992).
ity proportional to gradient. Although this algo-
Mathematically, the MFD algorithm can be
rithm produces fewer straight ¬‚ow pathways, it
written as
still restricts pathway segments to be multiples
of 45—¦ . Also, the fact that this algorithm does not p
max(0, Si )
fi = (3.1)
produce a unique answer is a problem for many j=1 (max(0, S j ))
In addition to the restricted orientations of where fi is the fraction of incoming ¬‚ow di-
the D8 and ρ8 algorithms, a larger problem ex- rected to each of the neighboring pixels, S is the
ists. Both of these models implicitly assume con- slope or gradient between the central point and
vergent ¬‚ow everywhere on the landscape. In each of its neighbors (indexed by i and j) and
nature, divergent ¬‚ow is common on hillslopes p is a free parameter. To calibrate this method,
and in distributary channel environments (e.g. Freeman (1991) used Eq. (3.1) to predict ¬‚ow on
alluvial fans) where ¬‚ow spreads out with dis- a conical surface for various values of p. For
p = 1.0, Freeman found that ¬‚ow was preferen-
tance downstream rather than becoming more
concentrated downstream. In the early days of tially directed towards diagonal pixels. Using a
slightly higher value of p = 1.1, this effect was
DEM analysis, many DEMs were so coarse (e.g.
eliminated. As such, p = 1.1 should be used for
90 m/pixel) that they did not resolve hillslopes or
channel beds. Single-direction ¬‚ow-routing algo- best results. Codes for implementing the MFD al-
rithms worked well on many of these early DEMs gorithm and for hydrologically correcting DEMs

(i.e. removing spurious pits and ¬‚ats) prior to ¬‚ow
routing are available in Appendix 2.
The difference between single and multiple-
direction ¬‚ow-routing algorithms is illustrated
in Figure 3.1. Figure 3.1b and 3.1c illustrate the
contributing area for Hanaupah Canyon, Death steepest descent MFD
Valley, California using the steepest-descent (D8) 117.00 116.95 116.90
and MFD methods, respectively. These algorithms 2 km
result in very different spatial distributions of 36.225
¬‚ow on the alluvial fan and on hillslopes, where
divergent ¬‚ow is important. Both ¬‚ow-routing 36.20
methods are approximations to the true distri-
bution of ¬‚ow, but the MFD method is more re-
alistic because it enables the ¬‚ow distribution to
depend continuously on the topographic shape steepest descent

while the steepest-descent method only depends (c)
on the drainage network topology. For example,
the steepest-descent method does not distinguish
between gentle, broad valleys and steep, narrow
valleys as long as they have the same drainage-
network topology. The MFD method, in contrast,
is sensitive to both topology and hypsometry. Not
surprisingly, the behavior of landform evolution MFD
models in which erosion and sediment transport
are a function of drainage area depends on the 106 107 m2
100 104
¬‚ow-routing algorithm used in the model (e.g.
Fig 3.1 Comparison of steepest-descent and MFD
Pelletier, 2004d).
¬‚ow-routing algorithms. (a) Steepest descent: a unit of
In addition to MFD, two additional multiple-
precipitation that falls on a grid point (shown at top) is
direction algorithms are widely used. Both of successively routed to the lowest of the eight nearest
these algorithms are similar to the steepest de- neighbors (including diagonals) until the outlet is reached.
scent algorithm in that they assume that ¬‚ow Precipitation can be dropped and routed in the landscape in
within a pixel has a well-de¬ned direction. How- any order. MFD: all incoming ¬‚ow to a grid point is
ever, rather than forcing the ¬‚ow direction to be distributed between the down-slope pixels, weighted by bed
a multiple of 45—¦ as in the steepest descent al- slope. To implement this algorithm ef¬ciently, grid points
should be ranked from highest to lowest, and routing should
gorithm, both the DEMON and D∞ algorithms
be done in that order to ensure that all incoming ¬‚ow from
compute a ¬‚ow direction that is allowed to ac-
up-slope is accumulated before downstream routing is
quire any value between 0 and 360—¦ . The D∞ al-
calculated. (b) Map of contributing area calculated with
gorithm (Tarboton, 1997) calculates the direction steepest descent for Hanaupah Canyon, Death Valley,
of steepest descent based on the slopes of eight California, using USGS 30 m DEMs. Grayscale is logarithmic
triangular facets determined by the central pixel and follows the legend at ¬gure bottom. (c) Map of
and pairs of adjacent neighboring pixels. Flow is contributing area computed with bifurcation routing, for the
then partitioned between two of the eight direc- same area and grayscale as (b). Multiple-¬‚ow-direction
routing results in substantially different and more realistic
tions nearest to the computed aspect, weighted
¬‚ow distribution, particularly for hillslopes and areas of
by their angular difference from the aspect. For
distributary ¬‚ow. For color version, see plate section.
example, if the ¬‚ow direction for a certain pixel is
Modi¬ed from Pelletier (2004d).
determined to be 10—¦ north of east, then 35/45ths
of the ¬‚ow is delivered to the pixel to the east
and 10/45ths of the ¬‚ow is delivered to the
pixel to the northeast. In the DEMON algorithm

creates too much ˜˜dispersion™™ (i.e. ¬‚ow spreads
0.8 km

out more quickly with downstream distance than
is realistic). There is no clear basis for conclud-
ing that MFD overpredicts the spread of ¬‚ow,
however. In creating his MFD algorithm, Freeman
(1991) began with the hypothesis that ¬‚ow on a
simple, symmetric landform such as a cone must
obey the symmetry of the cone. By assuming that

there is one principal direction of ¬‚ow for all pix-
els, the DEMON and D∞ algorithms violate that
symmetrical ¬‚ow hypothesis. On divides, for ex-
ample, where ¬‚ow tends to be directed roughly
equally in two opposite directions, the DEMON
and D∞ algorithms unrealistically force the ¬‚ow
to occur in one principal direction. The same bias
occurs in channels that bifurcate in two direc-
tions separated by more than 45—¦ .
Figure 3.3 compares the D∞ and MFD algo-
rithms for a 1 m/pixel DEM of a crater wall on
Mars. This landform is dominated by a steep
(20--35—¦ ) northwest-oriented slope dissected by
longitudinal channels 1 m deep. Both of these al-
gorithms yield similar spatial distributions for A,
but the MFD algorithm does create a somewhat
Fig 3.2 Maps of log(A ) estimated using the ¬ve algorithms
˜˜smoother™™ distribution, consistent with the fact
discussed in this chapter, using a high-resolution DEM of the
that ¬‚ow down such a slope will generally be
north ¬eld, an agricultural ¬eld in northeastern Colorado.
partitioned between three down-slope pixels (i.e.
From Erskine et al. (2006). Reproduced with permission of
the American Geophysical Union. north, northwest, and west) rather than only two
as in the D ∞ method. So, which map is cor-
rect? In order to answer that question, we should
(Costa-Cabral and Burges, 1994), grid elevation remember the goal of ¬‚ow-routing algorithms:
values are taken as pixel corners and ¬‚ow direc- to mimic the ¬‚ow of water (or some other con-
tions are based on the aspect of a plane surface stituent of mass or energy) through the land-
¬t to each pixel. As in the D∞ algorithm, ¬‚ow is scape. Viewed through this lens, we can con-
then partitioned into the two neighbors closest clude that both algorithms are limited by the
to the ¬‚ow direction. fact that they do not include any information
Several studies (e.g. Erskine et al., 2006) have on ¬‚ow depth. In relatively planar terrain such
conducted comparison tests of different ¬‚ow- as a steep crater slope, the spread of ¬‚ow across
routing algorithms used for computing con- and down the slope will depend on the ratio
tributing area. In general, these studies endorse of the ¬‚ow depth to the channel depth. If ¬‚ow
the use of multiple-¬‚ow-direction algorithms depths are comparable to or larger than chan-
over single-¬‚ow-direction algorithms for all appli- nel depths, ¬‚ow will be broadly distributed across
cations. Figure 3.2 illustrates that multiple-¬‚ow- the slope. Conversely, if ¬‚ow depths are small rel-
direction algorithms generate more realistic ¬‚ow ative to channel depths, ¬‚ow will be relatively
pathways than single-direction methods. concentrated in channels, as predicted by both
There is no widespread agreement regarding Figure 3.3b and 3.3c. In Section 3.4, we will dis-
which of the multiple-¬‚ow-direction algorithms cuss an algorithm based on the MFD model that
is best. In introducing his D∞ algorithm, Tar- successively ˜˜¬lls™™ the topography up to a speci-
boton (1997) argued that the MFD algorithm ¬ed ¬‚ow depth in order to quantify the effects of


shaded relief

(b) 200 m


low-relief terrain. Figure 3.4, for example, illus-
Fig 3.3 Grayscale maps of log(A ) estimated using (b) the
D∞ and (c) MFD algorithms, using a high-resolution DEM of trates a shaded-relief image of Hanaupah Canyon
a crater wall on Mars. Shaded relief image shown in (a), with fan in Death Valley. On the alluvial fan, where the
illumination from the northeast.
along-dip topography is locally planar, the topog-
raphy appears to have ˜˜bands™™ or ˜˜steps™™ in this
a ¬nite ¬‚ow depth on the spatial distribution of image. This banding is an artifact of the way that
¬‚ow. US Geological Survey DEMs are constructed.
USGS DEMs are made by interpolating 7.5 min
USGS contour maps. In high-relief areas, con-
3.3 “Cleaning up” US Geological tours are closely spaced and DEM values are
Survey DEMs quite accurate at the scale of 10 or 30 m. On
alluvial fans and other low-relief areas where
slopes are a few percent or less, contours can
The most widely used DEMs in the United States
be spaced by 100 m or more. Interpolating con-
are the US Geological Survey™s 10 m/pixel and
tours in areas with such low-density information
30 m/pixel DEMs. Despite their widespread use,
is problematic. DEMs have two kinds of problems
these DEMs suffer from well-known artifacts in

to the published contour map of the area. Then,
we encoded each of the contour pixels in the
digital contour map with elevation values from
contour artifacts

the DEM. The result is the elevation-coded raster
dataset of contours shown in Figure 3.6a.
The idea of the new interpolation method
is that variations in slope occur most gradually
along ¬‚ow lines. As such, the best estimate of
the elevation values between contours can be
obtained with a weighted interpolation along
¬‚ow lines. In order to perform this interpola-
tion, the distance between each pixel and its up-
stream and downstream contours must be com-
puted. Figure 3.6b presents a grayscale map of
the minimum distance to the downstream con-
tour. In general, ¬‚uvial topography evolves in or-
Fig 3.4 Shaded relief image of a portion of the USGS
der to minimize downstream transport distances.
30 m/pixel DEM of Hanaupah Canyon fan, Death Valley,
As such, ¬‚ow lines will generally correspond to
California, showing contour artifacts common in low-relief
pathways with the shortest distance down-slope.
To produce this minimum-downstream-distance
map, a local search algorithm was used that
calculates the distance between each pixel and
in low-relief areas. First, microtopographic vari-
all possible ¬‚ow pathways down-slope from that
ations will not be resolved in the original con-
pixel to the ¬rst contour line. The distance along
tour map, and hence they will also not be re-
each pathway was computed but only the mini-
solved in any DEM constructed from the contour
mum value was retained. In the up-slope direc-
map. Second, USGS DEMs are constructed by in-
tion it is more appropriate to use the maximum
terpolating between contours along N--S and E--W
distance to the upstream contour. This distance
grid lines. This interpolation method is simple
was calculated using a similar search algorithm
to implement but results in most of the banding
and is illustrated in Figure 3.6c.
present in USGS DEMs. A more natural interpola-
Figure 3.5b illustrates the DEM resulting from
tion method would interpolate between contours
a weighted linear interpolation of elevation val-
along ¬‚ow lines rather than along arti¬cial N--S
ues along ¬‚ow lines. The algorithm minimizes
and E--W grid lines. In this section, we explore the
some but not all of the contour artifacts present
use of ¬‚ow-line interpolation methods for ˜˜clean-
in the USGS DEM and generally provides for
ing up™™ some of the contour artifacts present in
crisper topographic detail. Figures 3.5a and 3.5c
can be compared to the ˜˜actual™™ topography in
Our study site is the Marshall Gulch water-
Figure 3.5c. This shaded-relief map was produced
shed near Summerhaven, Arizona. Figure 3.5a
from a 1 m/pixel LIDAR DEM of Marshall Gulch.
shows a shaded-relief image of the 10-m resolu-
Most of the small-scale detail present in this im-
tion USGS DEM of the Marshall Gulch watershed.
age is missing in both of the contour-derived
Contour artifacts are clearly present in this DEM.
DEMs. This underscores the fundamental limita-
In order to develop and test a new method for
tions of any DEM derived from relatively coarse
contour interpolation, we ¬rst need to backtrack
(20 ft) contour data. Clearly, LIDAR data is supe-
from the USGS DEM to the original contour data
rior where it is available. One reason why some
that was used to create it. We can do this by ¬rst
contour artifacts remain in the DEM is that a
creating a contour map from the DEM with the
linear interpolation along ¬‚ow lines, while supe-
same contour spacing as the original topographic
rior to a linear interpolation along N--S and E--W
map. The USGS contour map for this area has a
grid lines, nevertheless produces a slope break
20 ft contour interval. Using this interval, we cre-
at contour lines that may not be present in the
ated a digital contour map that is nearly identical

300 m

10 m/pixel USGS DEM 2 m/pixel flowline-interpolated DEM

(a) (b)

1 m/pixel LIDAR DEM
m/pixel LIDAR DEM
of water and water-borne constituents depend
Fig 3.5 Shaded-relief maps of Marshall Gulch DEMs.
on how much water is ¬‚owing across the land-
(a) USGS 10 m/pixel DEM, (b) ¬‚ow-line-interpolated
2 m/pixel DEM using method described in this section, scape. When the ¬‚ow depth in a channel ex-
(c) 1 m/pixel LIDAR DEM. ceeds bankfull capacity, for example, ¬‚ow will
be diverted onto the ¬‚oodplain. Since the MFD
and other multiple-direction ¬‚ow-routing meth-
actual topography. As such, the DEM of Figure
ods route ¬‚ow according to bed slope and do
3.5b could be further improved by using a spline
not incorporate ¬‚ow depth, they implicitly pre-
interpolation involving two or more contours in
dict ¬‚ow pathways for cases of negligible ¬‚ow
the upstream or downstream directions.
depth. In many geomorphic applications it is nec-
essary to constrain water and sediment pathways
3.4 Application of ¬‚ow-routing for a range of event sizes, not just for those cor-
responding to low ¬‚ow conditions. Furthermore,
algorithms to estimate ¬‚ood
in many cases modeling the detailed behavior of
hazards individual ¬‚ood events would be unrealistic. For
such applications it would be useful to have a
¬‚ow-routing method that retains the simplicity
A principal drawback of all of the ¬‚ow routing
of the MFD method yet is able to resolve the
methods we have considered is that they do not
different ¬‚ow pathways that occur when ¬‚ow
incorporate ¬‚ow depth. Clearly, the ¬‚ow paths

300 m

elevation-coded contour raster contour distance down flow line

(a) (b)

contour distance up flow line
Fig 3.6 Intermediate products used to create the MFD model is run again on the new partially
¬‚ow-line-interpolated DEM. (a) Elevation-encoded contour ¬lled DEM. This process is repeated using small
raster dataset, (b) grayscale map of minimum distance to
increments of ¬‚ow until the prescribed maxi-
downstream contour, (c) grayscale map of maximum distance
mum ¬‚ow depth is reached. This algorithm ef-
to upstream contour. Datasets in (b) and (c) provide
fectively mimics the ¬‚ood behavior of a slowly
weighting factors for interpolating between contours.
rising hydrograph. Done properly, this method
can be shown to converge to a unique solution
depths exceed the local cross-sectional area of the that satis¬es the ¬‚ow continuity equation as the
channel. number of iterations is increased.
The MFD model can be implemented in an it- Figure 3.7 illustrates the results of succes-
erative fashion to account for the effects of ¬‚ow sive MFD ¬‚ood routing on a high-resolution
depths on ¬‚ow pathways. The basic idea is to ¬rst (1 m/pixel) photogrammetric DEM of Fortymile
determine ¬‚ow pathways in low-¬‚ow conditions Wash in Amargosa Valley, Nevada. The initial con-
using the basic MFD routing method. Then, the dition for each of these runs is an upstream dis-
¬‚ow is rerouted on a DEM that includes the orig- charge prescribed at the location of the main
inal bed topography plus a small increment of channel at the northern boundary of the DEM.
¬‚ow determined by a prescribed maximum ¬‚ow The discharge value is ¬rst translated into a
depth together with the spatial distribution of ¬‚ow depth using Manning™s equation. Then, the
¬‚ow computed by the ¬rst MFD iteration. Then, MFD algorithm is run successively with larger

(a) (b)
flow depth


Q = 300 m3/s
1 km

(c) (d)

Q = 3000 m3/s
Q = 1000 m3/s

Fig 3.7 Illustration of a successive ¬‚ow-routing algorithm inundation, and can be used on a reconnoitering
applied to Fortymile Wash alluvial fan in Amargosa Valley, basis to determine relative ¬‚ood risk. Successive
Nevada. (a) Shaded-relief image of a high resolution
¬‚ow-routing analysis can be performed for a sin-
(1 m/pixel) photogrammetric DEM of Fortymile Wash, (b)“(d)
gle master channel network by specifying ¬‚ow
¬‚ow depths predicted by the model for different values of the
at an upstream boundary, or a uniform runoff
prescribed upstream discharge: (b) 300 m3 /s, (c) 1000 m3 /s,
depth can be speci¬ed for the entire DEM in or-
and (d) 3000 m3 /s. For color version, see plate section.
der to model distributed ¬‚ow routing. Appendix
2 provides code for implementing the distributed
increments of ¬‚ow until the input ¬‚ow depth ¬‚ow routing case.
corresponding to the upstream discharge is
achieved. This algorithm is useful at identifying
distributary channels and adjacent overbank ar-
eas that are activated when certain threshold 3.5 Contaminant transport in
¬‚ow depths are exceeded. This technique makes
channel bed sediments
a number of simplifying assumptions about the
steadiness of the ¬‚ow, and hence cannot be
used for detailed ¬‚ood-hazard assessment. Nev- Predicting the transport and fate of contami-
ertheless, it is a simple, fast method for iden- nants in ¬‚uvial systems is a critical aspect of
tifying ¬‚ow pathways and their thresholds of applied geomorphology. Worldwide, many river

systems have channel-bed sediments with in order to map the fraction of contaminated
elevated heavy-metal concentrations as a result sediments in all the channels of the Fortymile
of centuries of ore-mining activity (James, 1989; Wash drainage basin. The contaminant fraction
Miller, 1997). Similarly, radionuclide contami- computed for channels in the area of the critical
nation of channel-bed sediments is a potential population provides a basis for estimating the po-
human health hazard downstream of areas tential dose rate to that population in the event
of high-level nuclear processing and testing of an eruption.
(Reneau et al., 2004). The model divides the Fortymile Wash
In this section, we explore a model for the drainage basin into two domains: (1) the tribu-
prediction of contaminant distributions in ¬‚u- tary drainage basin of Fortymile Wash and (2) its
vial channel sediments. The motivation for de- distributary alluvial fan where the nearest pop-
veloping the model was the need to quantify ulation resides (Figure 3.8). The Fortymile Wash
the dispersal of radionuclide-contaminated sed- drainage basin includes the proposed repository
iments in the event of a volcanic eruption at location and the area of most of the primary fall-
the proposed nuclear-waste repository at Yucca out in the event of a volcanic eruption through
Mountain. Underground storage of nuclear waste the repository. The drainage basin and the fan are
primarily poses a hydrologic risk associated with divided at the fan apex. The model assumes that
contaminated groundwater. However, volcanism primary fallout is mobilized and transported to-
in the Yucca Mountain region could pose a signif- ward the alluvial fan if it falls on slopes steeper
icant geomorphic hazard if an eruption were to than a threshold gradient or on channels greater
occur. Probabilistic volcanic hazard analyses con- than a threshold stream power. To do this calcu-
strain the annual risk of an eruption through the lation, the model performs a spatially distributed
proposed repository to be a very low 1.5 — 10’8 analysis of the slopes and stream powers for the
(CRWMS M&O, 1996). In the event of an erup- entire Fortymile Wash drainage basin using an
tion, individuals living on the Fortymile Wash input 30-m resolution US Geological Survey Digi-
alluvial fan 18 km south of the repository (i.e. tal Elevation Model (DEM).
the closest residents to the repository; Figure 3.8) Before the mobilized tephra and radionu-
could be affected by radionuclide-contaminated clides are deposited on the alluvial fan, they are
tephra deposited as fallout or redistributed from transported through the alluvial channel system
upstream. Under most wind direction scenarios of Fortymile Wash where mixing with uncon-
(i.e. southerly winds), primary fallout tephra is taminated channel sediments results in the di-
concentrated to the north of the proposed repos- lution of contaminated tephra. During extreme
itory location and entirely within the Nevada Test ¬‚ood events, contaminated tephra deposited in
Site (BSC, 2004a,b). In such cases, ¬‚uvial redis- channels or mobilized into channels from ad-
tribution of contaminated tephra from the pri- jacent hillslopes will be entrained and poten-
mary fallout location to the Fortymile Wash allu- tially mixed with uncontaminated channel-bed
vial fan could be signi¬cant and hence must be sediments within the scour zone. The cumu-
evaluated. lative effect of many ¬‚oods is to mix the
In collaboration with scientists at Los Alamos tephra with any uncontaminated channel-bed
National Laboratory, my graduate students sediments located beneath the tephra but within
Stephen DeLong, Mike Cline, and I developed the scour zone both locally and with dis-
a model that quanti¬es the concentration of tance downstream. Quantifying this dilution pro-
tephra and radionuclides in channels of the cess requires a spatially distributed model of
Fortymile Wash alluvial fan as a result of hill- long-term ¬‚ood scour depths throughout the
slope and ¬‚uvial redistribution processes in the drainage basin. The MFD ¬‚ow routing model,
event of a volcanic eruption through the repos- together with empirical relationships between
itory. The model uses the MFD ¬‚ow-routing al- discharge and scour depth, provide the means
gorithm to route contaminated and uncontam- for mapping the scour depths necessary for the
inated sediments through the ¬‚uvial system model.

116.40 116.30
116.50° W

N 2 km
36.90 N
116° W
117° W
contours for
1 cm
southwesterly wind
1 mm
Nevada Fortymile
N Ca

ev lif
ad or

a ni

36.5 N Mtns.
Fortymile Wash

drainage boundary


N 10 km
36.0 N


Fig 3.8 Location map of the study areas in Amargosa Valley, non-source areas weighted by their relative up-
Nevada. Shaded relief image includes Lathrop Wells volcanic
stream areas:
center, the footprint of the proposed nuclear-waste
As A ns
C= Cs +
repository at Yucca Mountain, and the southern portion of C ns (3.2)
A s + A ns A s + A ns
the Fortymile Wash drainage basin (i.e. the site of primary
fallout for tephra with potential for redistribution to the where C is the contaminant concentration at a
Fortymile Wash fan where the nearest population resides. point downstream, A s and A ns are the source and
From Pelletier et al. (2008). Reproduced with permission of
non-source areas upstream from the point, and
Elsevier Limited.
C s and C ns are the source and non-source concen-
trations. The value of C ns is usually zero, but in
cases where contaminants have a ¬nite natural
Previous work on the dilution and mix- background level, or in cases where non-¬‚uvial
ing of ¬‚uvial-system contaminants predict one- (e.g. eolian) processes have transported contam-
dimensional (1D) concentration of contaminants inants from source to non-source areas in the
in channel systems. These models are 1D in basin, C ns may be greater than zero. Equation
the sense that vertical and cross-channel concen- (3.2) reduces to an exponential form if the source
tration variations are not considered. The clas- area is located in the basin headwaters and the
sic dilution-mixing model of Hawkes (1976) and upstream area is written as a function of x, the
Marcus (1987) assumes that the primary process along-channel distance:
of downstream dilution is the mixing of con-
C = C ns + C s exp(’x/xs ) (3.3)
taminants with uncontaminated sediments deliv-
ered from upstream. Mathematically, the model where xs is a length scale controlled by the
states that the concentration at a point down- relative areas of source and non-source regions
stream from a contaminant source is an average and the along-channel distance through the
of the contaminant concentration at source and source area. Equation (3.3) is an approximation

of Eq. (3.2); it describes the general downstream- con¬‚uence to be an average of the upstream
dilution trend but does not capture the abrupt concentrations, weighted by the upstream sedi-
concentration decrease that occurs as large trib- ment yields or basin areas. The model makes no
utaries enter the main channel. Equation (3.2) prediction regarding the vertical distribution of
has been tested in many contaminated rivers (e.g. contaminants.
Hawkes, 1976; Marcus, 1986) and its accuracy is The basic operation of the scour-dilution-
often remarkable considering the simplicity of mixing model is illustrated by the diagram in
the model relative to the complexity of the pro- the lower left corner of Figure 3.9b. The scour-
cesses involved. dilution-mixing model is conceptually similar to
In addition to mixing of the contaminant the classic model because the contaminant con-
with uncontaminated sediments from upstream, centration downstream is an average of upstream
dilution can also occur by mixing of contami- concentrations. However, the vertical component
nants with local channel-bed sediments down- of mixing is explicitly included, and the local
stream. In this process, a ¬‚ood wave moves contributions of contaminant and uncontami-
though the channel causing net erosion dur- nated channel-bed sediments are incorporated
ing the growing phase of the ¬‚ood, turbulent into the mixture. The basic assumption of the
mixing of the contaminant with uncontami- model is that the concentration at each point is
nated channel-bed sediments during the ¬‚ood, a mixture of contaminant and channel-bed sed-
and redeposition of the sediment--contaminant iments upstream from the point. As such, the
mixture as the ¬‚ood wave passes. The depth contaminant concentration at a channel point
of contaminant mixing by this process is the is a ratio of the total contaminant volume de-
scour depth, i.e. the thickness of sediments that livered from upstream to the total mobile-bed
undergo active transport during extreme ¬‚ood volume. One signi¬cant advantage of the scour-
events. The scour depth, in turn, is proportional dilution-mixing model is that it predicts the
to the square root of the unit discharge (dis- vertical distribution of contaminants as well as
charge per unit channel width) (Leopold et al., the surface concentration. In addition, by using
1966). A numerical model that maps unit dis- spatially distributed routing of runoff and con-
charge and uses Leopold et al.™s empirical rela- taminants, the model can predict cross-channel
tionship to predict scour depth, therefore, pro- distributions and can be applied to distributary
vides a basis for estimating the dilution effect environments. In this sense, the model predicts
caused by contaminant mixing with channel-bed a 3D distribution rather than the 1D distribution
sediments downstream. In our model, the mobile of the classic model.
portion of the bed (including both contaminants The scour-dilution-mixing model can be read-
and uncontaminated channel-bed sediments) is ily implemented using a raster-based framework.
assumed to be lifted from the bed, mixed by tur- In this framework, the contributing-area grid is
bulence, and deposited back onto the bed down- ¬rst initialized with the value of the pixel area,
( x)2 . Second, the MFD ¬‚ow-routing algorithm
stream. The dilution effect caused by this process
of Freeman (1991) is used to calculate the con-
can be quanti¬ed as a function of source and non-
tributing area routed through each pixel in the
source upstream areas as in the classic model.
DEM. Third, a critical stream-power threshold is
The classic dilution-mixing and scour-
used to distinguish between hillslopes and chan-
dilution-mixing models are compared schemat-
nels in the DEM, where the stream power is de-
ically in Figure 3.9. The basic operation of the
¬ned as the product of the local slope S and the
classic model is illustrated by the diagram
square root of the drainage area A (Montgomery
in the lower left corner of Figure 3.9a. The
and Dietrich, 1988). The critical stream power is
lengths of arrows coming into the central point
related to the drainage density X through the re-
represent the sediment yield from the two
lationship 1/ X . The value of the drainage density
upstream tributaries (one contaminated and
can be directly measured in the ¬eld or by using
one uncontaminated). The model predicts the
digital map products. Pixels with stream-power
contaminant concentration downstream of the

contaminant concentration




. 4
( 14)