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