<< . .

. 13
( : 51)

. . >>

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

<< . .

. 13
( : 51)

. . >>