<< . .

. 14
( : 51)



. . >>

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

<< . .

. 14
( : 51)



. . >>