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

2.3 NUMERICAL TECHNIQUES AND APPLICATIONS 57

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

3/2

Q x

c(x, 0, 0) = √ 1.0 m/kyr (assuming a Qa3 surface age of 50 ka),

(2.70)

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

area.

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-

applications

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-

method

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

expansion:

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

‚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)

58 THE DIFFUSION EQUATION

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

0

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

x of h or its ¬rst derivative must also be speci¬ed

x,t

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:

step.

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

‚x

2x x,t

a ¬nite value. The right-side boundary condition

(2.74)

is h n = 0. During the ¬rst time step, all of the

L

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

1

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

t

dreds of thousands of grid points, one million

Dtn

= 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

2.3 NUMERICAL TECHNIQUES AND APPLICATIONS 59

(b)

(a)

h h

unstable, ”t too large

h0

h0

h1

h1

stable

x

x

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

gives

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

0

=κ x 2 + (2.77)

explicit scheme, only the ¬rst grid point to the right, h 1 will

‚t ‚x ‚x

1

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

1

than h 1 . In the implicit scheme, the effect of raising the left sion of instantaneous base-level drop of an ini-

0

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.

60 THE DIFFUSION EQUATION

1.0

kt =

1.0

20 m2

0.8

0.8

kt = 0.6

0.6 h

h

1280 m2 h0

h0 kt =

1280 m2

0.4

0.4

kt =

kt =

q= k 0.2

0.2 2560 m2

2560 m2 kx

q= L

0.0

0.0

25

0 50

25

0 50

x (m)

x (m)

(a) (b)

kx

slope wash

1.0

q= L

without gulying

l

h

q= k

soil creep

h0

0.5

slope wash

with gullying

0.0

0.5

0 1.0

x/L

(c)

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

‚t

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

(2.78)

n

|‚h/‚ x|

1’ evolution decreases rapidly with the overall relief

Sc

2.3 NUMERICAL TECHNIQUES AND APPLICATIONS 61

evolution is very different from that of the linear

1.0

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

0.6

h broad range of ¬‚ux equations, initial conditions,

h0 and boundary conditions.

0.4

kt = 103 m2

2.3.4 2D Evolution of alluvial-fan terraces

0.2

Figure 2.23 illustrates solutions to the 2D diffu-

sion equation using the FTCS method. The model

0.0

30

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

(b)

(a)

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

remains.

t = 500 kyr

t = 120 kyr

62 THE DIFFUSION EQUATION

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-

κt

hi = hi + ’ 2h i + h i’1

n+1 n+1 n+1 n+1

n

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

⎜ ⎟

0

⎜ ⎟

’± (1 + 2±) ’± ...⎟

⎜0 no matter how large the grid is. This difference

⎜ ⎟

⎝0 ⎠

’± (1 + 2±) makes the implicit method far more stable than

0

... 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.

2

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

κt

h in+1 = h in + h i+1 ’ 2h in+1 + h i’1

n+1 n+1

’1

=A h

n+1 n

h (2.83) 2( x)2

+ h i+1 ’ 2h in + h i’1

n n

(2.84)

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.

EXERCISES 63

(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-

method

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-

κt

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

n

2

( 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-

(2.85)

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

Exercises

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

n+1

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

64 THE DIFFUSION EQUATION

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 .

(b)

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

EXERCISES 65

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-

3.2 ALGORITHMS 67

(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)

p

8

produce a unique answer is a problem for many j=1 (max(0, S j ))

applications.

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

68 FLOW ROUTING

(i.e. removing spurious pits and ¬‚ats) prior to ¬‚ow

(a)

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

(b)

and MFD methods, respectively. These algorithms 2 km

result in very different spatial distributions of 36.225

N

¬‚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-

36.175

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

102

¬‚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

3.2 ALGORITHMS 69

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

r8

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

70 FLOW ROUTING

(a)

shaded relief

(c)

(b) 200 m

MFD

D-infinity

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

3.3 “CLEANING UP” US GEOLOGICAL SURVEY DEMS 71

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.

areas.

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

USGS DEMs.

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

72 FLOW ROUTING

300 m

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

(a) (b)

1 m/pixel LIDAR DEM

m/pixel LIDAR DEM

(c)

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

3.4 APPLICATION OF FLOW-ROUTING ALGORITHMS TO ESTIMATE FLOODHAZARDS 73

300 m

elevation-coded contour raster contour distance down flow line

(a) (b)

contour distance up flow line

(c)

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

74 FLOW ROUTING

(a) (b)

flow depth

1m

0.3

0.1

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

3.5 CONTAMINANT TRANSPORT IN CHANNEL BED SEDIMENTS 75

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.

76 FLOW ROUTING

116.40 116.30

116.50° W

N 2 km

36.90 N

hypothetical

116° W

117° W

tephra-thickness

Yucca

contours for

1 cm

Mtn.

southwesterly wind

1 mm

repository

Nevada Fortymile

footprint

Wash

N Ca

Amargosa

ev lif

ad or

Valley

Funeral

36.80

a ni

36.5 N Mtns.

Fortymile Wash

a

drainage boundary

Death

Valley

Black

Mtns.

N 10 km

36.0 N

36.70

FortymileWash

fan

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

3.5 CONTAMINANT TRANSPORT IN CHANNEL BED SEDIMENTS 77

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

78 FLOW ROUTING

contaminant concentration

(a)

(b)

+

H