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