A new method for the determination of ﬂow directions
and upslope areas in grid digital elevation models
David G. Tarboton
Utah Water Research Laboratory, Utah State University, Logan
Abstract. A new procedure for the representation of ﬂow directions and calculation of
upslope areas using rectangular grid digital elevation models is presented. The procedure
is based on representing ﬂow direction as a single angle taken as the steepest downward
slope on the eight triangular facets centered at each grid point. Upslope area is then
calculated by proportioning ﬂow between two downslope pixels according to how close this
ﬂow direction is to the direct angle to the downslope pixel. This procedure offers
improvements over prior procedures that have restricted ﬂow to eight possible directions
(introducing grid bias) or proportioned ﬂow according to slope (introducing unrealistic
dispersion). The new procedure is more robust than prior procedures based on ﬁtting
local planes while retaining a simple grid based structure. Detailed algorithms are
presented and results are demonstrated through test examples and application to digital
elevation data sets.
Introduction
Flow directions based on digital elevation models (DEMs)
are needed in hydrology to determine the paths of water,
sediment, and contaminant movement. Two important distributed
quantities that depend on ﬂow directions are the upslope
area and speciﬁc catchment area. Upslope area, A, is deﬁned
as the total catchment area above a point or short length of
contour [Moore et al., 1991]. The speciﬁc catchment area, a, is
deﬁned as the upslope area per unit width of contour, L, (a ϭ
A/L) [Moore et al., 1991] and is a distributed quantity that has
important hydrological, geomorphological, and geological signiﬁcance
[Costa-Cabral and Burges, 1994]. The speciﬁc catchment
area contributing to ﬂow at any particular location is
useful for determining relative saturation and generation of
runoff from saturation excess in models such as TOPMODEL
[Beven and Kirkby, 1979; Beven et al., 1984; Wood et al., 1990].
Speciﬁc catchment area together with other topographic parameters
has also been used in the analysis of processes such as
erosion and landslides [Dietrich et al., 1992, 1993; Wu, 1993;
Montgomery and Dietrich, 1994]. Upslope area is commonly
used for the automatic demarcation of channels relying on the
notion of a critical support area [O’Callaghan and Mark, 1984;
Jenson and Domingue, 1988; Morris and Heerdegen, 1988; Tarboton,
1989; Lammers and Band, 1990; Tarboton et al., 1991,
1992; Martz and Garbrecht, 1992]. Judging by the number of
recent papers, there is considerable hydrologic interest in the
effect of grid scale and procedures for computation of speciﬁc
catchment area [Zhang and Montgomery, 1994; Quinn et al.,
1995; Wolock and McCabe, 1995]. It is therefore important that
ﬂow directions and speciﬁc catchment areas be accurately determined
free from grid artifacts.
In this paper a new method for calculating ﬂow directions,
upslope, and speciﬁc catchment areas is presented. In what
follows I ﬁrst review the currently available procedures for
calculating upslope and speciﬁc catchment areas on the basis
of grid digital elevation models, pointing out their strengths
and weaknesses and giving the motivation for a new method. I
then describe the new procedure for calculation of ﬂow directions
and the procedure for calculation of upslope area on the
basis of the new ﬂow directions. Illustrative examples are used
to compare the new method with existing methods. These
include simple test cases, where the theoretical result is known
and bias and square error can be evaluated, as well as low- and
high-resolution DEM data where the evaluation is qualitative.
Readers familiar with the issues surrounding the computation
of upslope area based on grid DEMs could skip the background
section.
Background
Grid DEMs consist of a matrix data structure with the topographic
elevation of each pixel stored in a matrix node. Grid
DEMs are distinct from other DEM representations such as
triangular irregular network (TIN) and contour-based data
storage structures. Grid DEMs are readily available and simple
to use and hence have seen widespread application to the
analysis of hydrologic problems [Moore et al., 1991]. However,
they suffer from some drawbacks that arise from their grided
nature.
The earliest and simplest method for specifying ﬂow directions
is to assign ﬂow from each pixel to one of its eight
neighbors, either adjacent or diagonal, in the direction with
steepest downward slope. This method, designated D8 (eight
ﬂow directions), was introduced by O’Callaghan and Mark
[1984] and has been widely used [Marks et al., 1984; Band,
1986; Jenson and Domingue, 1988; Mark, 1988; Morris and
Heerdegen, 1988; Tarboton et al., 1988; Tarboton, 1989; Jenson,
1991; Martz and Garbrecht, 1992]. In the context of a grid the
upslope area, A, is the area contributing to each pixel and may
be estimated as the product of the number of pixels draining
through each pixel and pixel area. The speciﬁc catchment area,
a, is then estimated as A/L, taking L as the pixel width.
The D8 approach has disadvantages arising from the discretization
of ﬂow into only one of eight possible directions,
separated by 45Њ [e.g., Fairﬁeld and Leymarie, 1991; Quinn et
Copyright 1997 by the American Geophysical Union.
Paper number 96WR03137.
0043-1397/97/96WR-03137$09.00
WATER RESOURCES RESEARCH, VOL. 33, NO. 2, PAGES 309–319, FEBRUARY 1997
309
al., 1991; Costa-Cabral and Burges, 1994]. Fairﬁeld and Leymarie
[1991] suggested overcoming these problems by randomly
assigning a ﬂow direction to one of the downslope
neighbors, with the probability proportional to slope.
Multiple ﬂow direction methods [Quinn et al., 1991; Freeman,
1991] have also been suggested as an attempt to solve the
limitations of D8. These allocate ﬂow fractionally to each lower
neighbor in proportion to the slope (or in the case of Freeman’s
method, slope to an exponent) toward that neighbor.
Multiple ﬂow direction methods, designated here by MS (multiple
directions based on slope), have the disadvantage that
ﬂow from a pixel is dispersed to all neighboring pixels with
lower elevation.
Dispersion is inherent in any method (including the one I
describe below) that assigns ﬂow from one pixel to more than
one downslope neighbor and manifests itself in terms of
spreading of ﬂow from a single pixel. It could be argued that
this does not matter because the models that use a may use it
as a surrogate for a physical quantity that is affected by dispersion.
However, dispersion is inconsistent with the physical definition
of upslope area, A, and speciﬁc catchment area, a. It is
important, to the extent possible, to minimize dispersion in the
calculation of a. Then, if necessary, physical dispersion can be
modeled separately.
Lea [1992] developed an algorithm that uses the aspect associated
with each pixel to specify ﬂow directions. Flow is
routed as though it were a ball rolling on a plane released from
the center of each grid cell. A plane is ﬁt to the elevations of
pixel corners, these corner elevations being estimated by averaging
the elevations of adjoining pixel center elevations. This
procedure has the advantage of specifying ﬂow direction continuously
(as an angle between 0 and 2␲) and without dispersion.
Costa-Cabral and Burges [1994] presented an elaborate
set of procedures, named DEMON (digital elevation model
networks), that extends the ideas of Lea [1992]. Grid elevation
values are used as pixel corners, rather than block centered,
and a plane surface is ﬁtted for each pixel. Costa-Cabral and
Burges [1994] recognized ﬂow as two-dimensional originating
uniformly over the pixel area, rather than tracking ﬂow paths
from the center point of each pixel. They evaluate upslope area
through the construction of detailed ﬂow tubes.
The assumption of a plane ﬁt locally to each pixel requires
approximation because only three points are required to determine
a plane. The best ﬁt plane cannot in general pass
through the four corner elevations, leading to a discontinuous
representation of the surface at pixel edges. Planes ﬁt locally to
certain elevation combinations can lead to inconsistent or
counterintuitive ﬂow directions that are a problem in both
Lea’s [1992] method and in DEMON. Figure 1 illustrates some
of these problems. These problems illustrated in the context of
Lea’s method are also present in DEMON, since the same
plane ﬂow directions would arise given the corner elevations
shown in Figure 1. The coding of approaches based on ﬁtting
local planes, such as Lea’s method and DEMON, so that they
are robust and work for all possible elevation combinations
that may arise in real data is difﬁcult. There are many exceptions,
such as the one illustrated in Figure 1, that need to be
anticipated and special code developed to account for them. In
fact, the code for DEMON upslope [Costa-Cabral and Burges,
1994] is unavailable because it is “hard to program and full of
special cases” (M. Costa-Cabral, personal communication,
1995).
Figure 1. Hypothetical DEM subset on a block-centered grid. Corner elevations are calculated by averaging
the elevations of adjoining cells. (a) Flow directions determined using locally ﬁtted planes [Lea, 1992]. (b)
Flow directions by this papers new procedure. In Figure 1a consider ﬂow from the pixel with center marked
A. The pixel immediately to the north is the lowest neighbor, with elevation 5, so that is where water from this
pixel should ﬂow. Some water from pixel A should also ﬂow to the south or south west, since pixel A is a
saddle. Lea’s method has water ﬂow to the northwest along the direction of the plane slope to point E, where
it crosses into pixel B. The edge between pixels A and B is a discontinuous edge with ﬂow from the adjacent
pixels in opposing directions, a situation that poses problems for plane ﬂow algorithms. Flow from A is routed
south along the edge from E to F, where it crosses into pixel C and without exception checking gets stuck in
an inﬁnite loop, moving from pixel C to D to A to B around the point F. Furthermore the ﬂow originating at
the center of pixel B goes south east to point H, where it enters pixel A and then ﬂows north along the edge
(since the slope vector centered at A has a north component), to point I, in an uphill direction, where it enters
the pixel J. Thus ﬂows from B and A move in opposite directions along the same edge, an unrealistic situation.
Another problem is evident at pixel K. Here ﬂow has a large north component due to the elevation 30 to the
south, despite the fact that its lowest neighbor is to the southeast, and neighbor with steepest slope is to the
east.
TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS310
The motivation for developing a new method stemmed from
frustration with the way these existing procedures worked. The
following issues are relevant to the evaluation and design of
DEM ﬂow direction procedures: (1) the need to avoid or
minimize dispersion; (2) the need to avoid grid bias, due to the
orientation of the numerical grid; (3) the precision with which
ﬂow directions are resolved; (4) a simple and efﬁcient grid
based matrix storage structure; and (5) robustness, the ability
to cope with “difﬁcult” data, such as the saddle in Figure 1 but
also with pits and ﬂat areas.
The D8 method does well on points 1, 4, and 5 but resolves
ﬂow directions too coarsely (point 3), introducing grid bias
(point 2). The randomness in the method of Fairﬁeld and
Leymarie [1991] is not appealing. Upslope and speciﬁc catchment
areas are deterministic quantities that we should be able
to compute in a repeatable way. The MS method avoids grid
bias (point 2) but introduces substantial dispersion. Since ﬂow
may be proportioned in up to eight possible directions, eight
numbers need to be stored (or recalculated each time they are
needed) for each pixel, resulting in inefﬁcient data storage
(point 4). The plane ﬂow methods [Lea, 1992; Costa-Cabral
and Burges, 1994] are appealing in that they are deterministic
and precisely resolve ﬂow directions. However, they are susceptible
to problems (point 5) arising from the approximation
involved in ﬁtting a plane through four points.
Calculation of Flow Directions
Here I propose a new method for the representation and
calculation of ﬂow directions that is a compromise on the
issues raised in the previous section. It recognizes the advantages
of Lea’s [1992] method and DEMON through the assignment
of a single ﬂow direction to each cell. This single ﬂow
direction (represented as a continuous quantity between 0 and
2␲) is determined in the direction of the steepest downwards
slope on the eight triangular facets formed in a 3 ϫ 3 pixel
window centered on the pixel of interest. The use of triangular
facets avoids the approximation involved in ﬁtting a plane and
the inﬂuence of higher neighbors on downslope ﬂow. Where
the direction does not follow one of the cardinal (0, ␲/2, ␲, and
3␲/2) or diagonal (␲/4, 3␲/4, 5␲/4, and 7␲/4) directions, upslope
area is calculated by apportioning the ﬂow from a pixel
between the two downslope pixels according to how close the
ﬂow angle is to the direct angle to that pixel center. Since only
a single number need be saved for each pixel to represent the
ﬂow ﬁeld, computer storage is simple and efﬁcient. Some dispersion
is introduced by the proportioning of ﬂow between
downslope pixels, but this is minimized since ﬂow is never
proportioned between more than two downslope pixels.
Figure 2 illustrates the calculation of ﬂow directions. A
block-centered representation is used with each elevation
value taken to represent the elevation of the center of the
corresponding pixel. Eight planar triangular facets are formed
between the pixel and its eight neighbors. Each of these has a
downslope vector which when drawn outward from the center
may be at an angle that lies within or outside the 45Њ (␲/4
radian) angle range of the facet at the center point. If the slope
vector angle is within the facet angle, it represents the steepest
ﬂow direction on that facet. If the slope vector angle is outside
a facet, the steepest ﬂow direction associated with that facet is
taken along the steepest edge. The ﬂow direction associated
with the pixel is taken as the direction of the steepest downslope
vector from all eight facets.
To implement this procedure ﬁrst consider a single triangular
facet (Figure 3). Slope (downward) is represented by the
vector (s1, s2) where
s1 ϭ ͑e0 Ϫ e1͒/d1 (1)
s2 ϭ ͑e1 Ϫ e2͒/d2 (2)
Figure 2. Flow direction deﬁned as steepest downward slope
on planar triangular facets on a block-centered grid.
Figure 3. Deﬁnition of variables for the calculation of slope
on a single facet.
Table 1. Facet Elevation and Factors for Slope and Angle Calculations
Facet 1 2 3 4 5 6 7 8
e0 ei,j ei,j ei,j ei,j ei,j ei,j ei,j ei,j
e1 ei,jϩ1 eiϪ1,j eiϪ1,j ei,jϪ1 ei,jϪ1 eiϩ1,j eiϩ1,j ei,jϩ1
e2 eiϪ1,jϩ1 eiϪ1,jϩ1 eiϪ1,jϪ1 eiϪ1,jϪ1 eiϩ1,jϪ1 eiϩ1,jϪ1 eiϩ1,jϩ1 eiϩ1,jϩ1
ac 0 1 1 2 2 3 3 4
af 1 Ϫ1 1 Ϫ1 1 Ϫ1 1 Ϫ1
311TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS
and ei and di are elevations and distances between pixels as
labeled in Figure 3. The slope direction and magnitude are
r ϭ tanϪ1
͑s2/s1͒
(3)
s ϭ ͑s1
2
ϩ s2
2
͒1/ 2
If r is not in the range (0, tanϪ1
(d2/d1)), then r needs to be
set as the direction along the appropriate edge and s assigned
as the slope along that edge.
if r Ͻ 0, r ϭ 0, s ϭ s1 (4)
if r Ͼ tanϪ1
͑d2/d1͒, r ϭ tanϪ1
͑d2/d1͒,
(5)
s ϭ ͑e0 Ϫ e2͒/͑d1
2
ϩ d2
2
͒1/ 2
Next recognize that each of the eight facets depicted in Figure
2 can be mapped by appropriate selection of corner elevations
and rotation/transformation onto the facet in Figure 3. Table 1
gives the node elevations corresponding to the corners of each
of the triangular facets used to calculate slopes and angles in
(1)–(5). These are arranged such that e0 is the center point, e1
is the point to the side, and e2 is the diagonal point. The local
angle associated with the largest downwards slope from the
eight facets (rЈ ϭ r from facet with maximum s) is then
adjusted to reﬂect an angle counterclockwise from east (Figure
2) to obtain the ﬂow direction angle.
rg ϭ af rЈ ϩ ac␲/ 2 (6)
The multiplier af and constant ac depend on the facet selected
and are listed in Table 1. The procedure that searches for the
facet with the largest slope proceeds in the order of facets 1 to
8 shown in Figure 1 and in the case of ties (facets with equal
slope) picks the ﬁrst one. In nature ties are extremely rare so
the bias introduced by this is deemed negligible.
In the case where no slope vectors are positive (downslope)
a ﬂow direction angle of Ϫ1 is used to ﬂag the pixel as “unresolved,”
that is, a ﬂat area or pit. Unresolved ﬂow directions
are resolved iteratively by making them ﬂow toward a neighbor
of equal elevation that has a ﬂow direction resolved. This is the
same approach for resolving pits and ﬂats as used in the D8
method [e.g., Mark, 1988; Jenson and Domingue, 1988]. I therefore
use the calculation of D8 ﬂow directions as a preprocessor
to raise the elevation of all pixels in a pit to the level of the
overﬂow. Then where pixels are ﬂagged as “unresolved,” the
ﬂow angle returned by the D8 procedure is used. This ensures
that ﬂat pixels drain to a neighbor that ultimately drains to a
lower elevation, eliminating the possibility of inconsistencies
such as loops in the ﬂow direction angles. This method of representing
ﬂow directions based on triangular facets is designated
Dϱ (an inﬁnite number of possible single ﬂow directions).
Figure 4. Top half of an outward draining circular cone. Elevation was deﬁned as 200 minus the radius from
the center on a 16 ϫ 16 grid with grid spacing 10 units. Speciﬁc catchment area is theoretically radius/20
ranging from 0 at the center to 53 on the corners. Contours (10-unit interval) depict elevation. Gray scale
depicts contributing area (1 white to 60 black). (a) Theoretical values, (b) single direction (D8) procedure, (c)
Quinn et al.’s [1991] procedure (MS), (d) Lea’s [1992] method, (e) DEMON [Costa-Cabral and Burges, 1994],
and (f) new procedure (Dϱ).
TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS312
Calculation of Upslope Areas
Upslope area is calculated using a recursive procedure that
is an extension of the very efﬁcient recursive algorithm for
single directions [Mark, 1988]. The upslope area of each pixel
is taken as its own area (one) plus the area of upslope neighbors
that have some fraction draining to the pixel in question.
The ﬂow from each cell either all drains to one neighbor (if the
angle falls along a cardinal or diagonal direction) or is on an
angle falling between the direct angle to two adjacent neighbors.
In the latter case the ﬂow is proportioned between these
two neighbor pixels according to how close the ﬂow direction
angle is to the direct angle to those pixels. The following
pseudocode gives the logic of this algorithm:
Procedure DPAREA(i, j)
if AREA(i, j) is known
then
no action
else
AREA(i, j) ϭ 1 (the area of a single pixel)
for each neighbor (location in, jn)
p ϭ proportion of neighbor (in, jn) that drains to
pixel (i, j) based on angle
if ( p Ͼ 0) then
call DPAREA(in, jn) (this is the recursive call to
calculate area for the neighbor)
AREA(i, j) ϭ AREA(i, j) ϩ p ϫ
AREA(in, jn)
Return
The calculation is initiated by calling this function for the
outlet pixel. It then recursively calls itself for all pixels that
contribute to the upslope area at the outlet. The recursion
stops when it reaches a pixel that has no pixels upslope.
Illustrative Examples
This section gives examples of results from this method, Dϱ,
compared to the single direction approach, D8; Quinn et al.’s
[1991] multidirection algorithm, MS; Lea’s [1992] method; and
DEMON [Costa-Cabral and Burges, 1994]. In these examples
we use the notion of inﬂuence and dependence maps. The
inﬂuence function I( x, x0) is deﬁned as the upslope area at
each pixel x from a speciﬁc pixel x0. It maps where ﬂow from
pixel x0 goes and how it is dispersed. It is computed by running
a modiﬁed version of the procedure for calculating upslope
area that uses an area contribution of one from pixel x0 but
zero for all other pixels. The dependence function D( x, x0) is
the opposite of the inﬂuence function, deﬁned as D( x, x0) ϭ
I( x0, x). The upslope area at pixel x0 is composed of the sum
of the area of upslope pixels that have some proportion of their
ﬂow go through pixel x0. D( x, x0) maps the contribution from
pixel x to the calculation of upslope area at x0. It is calculated
through repeated evaluation of the inﬂuence function.
Figure 4 shows the upslope area by each approach for a
circular cone. Figure 5 shows the inﬂuence maps from each of
the ﬁve algorithms (D8, MS, Lea’s [1992] method, DEMON,
and Dϱ) applied to the circular cone. D8 results in no spreading,
but ﬂow paths (which are what the inﬂuence map plots)
Figure 5. Circular cone inﬂuence maps for the circled pixels. Gray scale ranges from white (0 or no
inﬂuence) to black (1 or 100% inﬂuence).
313TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS
are constrained to grid directions. MS follows the topographic
slope but introduces substantial dispersion. Lea’s [1992]
method marches down the contours in stair step fashion, while
DEMON and Dϱ strike a balance with spreading slightly wider
than divergence perpendicular to the contours. For some pixels
(e.g., the lower right one in Figure 5e) Dϱ results in no spreading.
This is a pixel where the ﬂow direction is aligned with the
grid diagonal. In general, when the topographic slope is
Figure 6. Inward cone. Contours show elevation. The left
panels show upslope area in gray scale. The right panels show
the inﬂuence function from the circled pixel. (a) D8 method,
(b) MS method [Quinn et al., 1991], (c) Lea’s [1992] method,
(d) DEMON [Costa-Cabral and Burges, 1994], and (e) Dϱ
method.
Figure 7. Dependence maps for planar surface. Gray scale (0
white, 1 black) indicates the fraction of each pixel upslope to
the circled pixel.
TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS314
Figure 8. Images of upslope area for a portion of the Walnut Gulch Experimental Watershed DEM. A
logarithmic scale of gray shades is used with lighter shades corresponding to higher values. This is a U.S.
Geological Survey DEM with 30-m resolution. (a) Single-direction, D8, procedure; (b) Quinn et al. [1991]
procedure, MS; and (c) new procedure, Dϱ.
315TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS
aligned with the grid axes, cardinal or diagonal, the Dϱ procedure
gives the same results as D8, and both are correct.
However, when the topographic slope is not aligned with one
of the grid directions, the procedures differ. D8 introduces no
dispersion, but at the expense of grid bias. Dϱ follows the
topographic slope at the cost of introducing some dispersion.
Figure 6 shows upslope area and the inﬂuence functions for
a portion of an inward ﬂowing cone surface. Figure 7 shows the
dependence maps for a plane surface not aligned with the grid.
The dependence map reﬂects the fraction of a pixel’s area that
drains to the designated pixel. It serves to demarcate a zone of
contribution, with shading to denote the degree of contribution.
On a planar surface the dependence maps should be
straight lines perpendicular to the gradient. The D8 method
gives straight lines following grid lines. MS has dependence
from a broad area, illustrating strikingly the problems with MS,
even for a simple surface. The dependence map for Lea’s
[1992] method is a stair step path roughly perpendicular to the
contours, as expected. DEMON is similar with some spillover
into adjacent pixels due to the two-dimensional ﬂow representation.
Dϱ has dependence from a narrower band 45Њ wide
upslope of the point under consideration. The axis of this band
is perpendicular to the contours.
For each of these example surfaces, cone, inward cone, and
plane, the true upslope area was computed at each grid point
and compared to results from each of the DEM procedures.
Table 2 presents the differences between the theoretical result
and each DEM procedure.
In evaluating the error statistics the grid direction bias introduced
by D8, clearly evident in the ﬁgures, is responsible for
the large mean square error (MSE) on the cone surfaces.
Curiously, D8 does well for the plane because the area is the
same whether one counts along the grid or perpendicular to
contours. This would not have been the case had the ridge not
been aligned with the grid. Quinn et al.’s [1991] MS method
does best for the inward cone where the concave surface limits
Figure 8. (continued)
Table 2. Differences Between Theoretical and DEM-Computed Upslope Area for Test Examples Expressed in Terms of
the Mean Error and Mean Square Error
Outward Cone Inward Cone Plane
Bias
Mean ( A Ϫ Aˆ )
MSE
Mean (( A Ϫ Aˆ )2
)
Bias
Mean ( A Ϫ Aˆ )
MSE
Mean (( A Ϫ Aˆ )2
)
Bias
Mean ( A Ϫ Aˆ )
MSE
Mean (( A Ϫ Aˆ )2
)
D8 Ϫ0.13 2.13 1.76 118.88 Ϫ0.17 0.065
MS Ϫ0.81 0.69 Ϫ1.07 5.70 Ϫ1.37 2.065
Lea’s [1992] method Ϫ1.29 2.41 Ϫ4.05 44.00 Ϫ2.57 7.912
DEMON Ϫ0.37 0.17 Ϫ0.37 19.23 Ϫ0.40 0.161
Dϱ Ϫ0.13 0.20 1.87 30.58 Ϫ0.17 0.065
Bias, mean error; MSE, mean square error; A, true upslope area; Aˆ, computed upslope area.
TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS316
dispersion. However, it does poorly for the plane and outward
cone owing to its substantial lateral dispersion. For example,
with the outward cone, MS pixels along the outward edge have
a component of slope parallel to the edge towards the corners.
Therefore some fraction of ﬂow that should exit the domain at
the center is directed towards the corners, increasing the sum
of upslope area around the edge above the area of the domain.
Lea’s [1992] method in all cases has large MSE and bias due to
a bias toward overestimating upslope area by counting a contribution
of 1 even if a ﬂow path from a pixel only just crosses
the corner of a pixel. In terms of statistics DEMON and Dϱ
are quite competitive. Both have small MSE and bias for the
outward cone and plane and small MSE relative to Lea’s and
D8 for the inward cone. DEMON is better on the inward cone.
This is due to the continuity across pixels maintained by DEMON
with its ﬂow tube procedure. Toward the middle of the
Figure 9. Images of upslope area for a portion of the Tennessee Valley study area DEM, Marin County,
California. A logarithmic scale of gray shades is used with lighter shades corresponding to higher values. This
is a 2-m resolution DEM generated from low-altitude stereo aerial photographs [Dietrich et al., 1992, 1993;
Montgomery and Foufoula-Georgiou, 1993]. (a) Single-direction, D8, procedure; (b) Quinn et al. [1991]
procedure, MS; and (c) new procedure, Dϱ.
317TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS
cone DEMON has many ﬂow tubes at spacing much ﬁner than
the cell size. Dϱ mixes the ﬂow in each cell resulting in some
blockiness evident in Figure 6e. Though informative, these
statistics do not provide the complete picture. The good performance
of MS on the inward cone is due to its dispersion
being contained and considerable concentric averaging
smoothing grid effects that are not so smoothed by the other
methods. However, this averaging which contributes to small
error statistics for a symmetric test case is in general undesirable
as it results in loss of resolution.
Figure 8 shows the upslope area for portion of the Walnut
Gulch Experimental Watershed calculated using D8, MS, and
Dϱ. This is a 30-m resolution DEM from the U.S. Geological
Survey Hay Mountain and Tombstone 7.5-min quadrangles.
The differences at this scale are quite subtle, but it is possible
to see more streakiness associated with the grid from the D8
procedure (Figure 8a) than with that from MS (Figure 8b) and
Dϱ (Figure 8c).
Figure 9 shows the upslope area for a portion of the Tennessee
Valley area near San Francisco, California. This is a
high-resolution (2-m grid) DEM generated from low-altitude
stereo aerial photographs [Dietrich et al., 1992, 1993; Montgomery
and Foufoula-Georgiou, 1993]. At this scale the differences
are much more noticeable than was the case for Figure 8. The
D8 procedure results in streaks aligned with the grid indicating
that ﬂow down hillsides is biased by the grid alignment. The
Dϱ and MS procedures do not suffer from these problems. The
Dϱ results seem to have more sharpness than MS. This is due
to less grid dispersion. Upslope areas were not computed with
DEMON or Lea’s [1992] method for the real DEMs, because
in the case of Lea’s method it performed poorly in the test
example evaluation (Table 2), and the DEMON code does not
work automatically given the ﬂat areas and situations requiring
exceptions present in this data. The DEMON downslope code
identiﬁed 6607 sinks in the Walnut Gulch DEM, even after the
pits had been ﬁlled.
In Figure 10, histograms showing the distribution of speciﬁc
catchment area by the methods D8, MS, and Dϱ are shown.
These indicate signiﬁcant differences in the number of pixels
computed to have a certain speciﬁc catchment area at the low
end of the speciﬁc catchment area scale. At the high end the
methods all give similar results. Speciﬁc catchment area is
Figure 9. (continued)
Figure 10. Histograms indicating the distribution of speciﬁc
catchment area by each method. (a) Portion of Walnut Gulch
Watershed shown in Figure 8. (b) Portion of Tennessee Valley
area shown in Figure 9.
TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS318
small on hillslopes and large in valleys and on the channel
networks. Therefore, if the intent is demarcation of valleys or
channels, for example, using the notion of a critical support
area, the different methods will give similar results. However,
if the intent is to use the speciﬁc catchment areas for distributed
modeling of hydrologic processes, such as runoff generation
or erosion, then the different methods will give differing results.
Conclusions
A new procedure for the representation of ﬂow directions
and the calculation of upslope area using grid-based digital
elevation models was presented. The procedure is based on
representing ﬂow direction as a single angle taken as the steepest
downward slope on the eight triangular facets centered at
each pixel. Results from the procedure were compared to
other grid-based procedures for calculation of upslope area
from grid DEMs, using test examples and low- and highresolution
DEM data. On the basis of the evaluation of test
statistics and examination of inﬂuence and dependence maps,
the new procedure performs better than D8, Lea’s [1992]
method, and MS and is comparable to DEMON [Costa-Cabral
and Burges, 1994]. The new procedure overcomes the problems
of loops and inconsistencies that plague plane-ﬁtting methods
such as DEMON. In real DEMs, signiﬁcant differences between
the distribution of speciﬁc catchment area were obtained
depending on the method used. Differences are more
marked as digital elevation data resolution increases, especially
at hillslope scales. Where results from the different methods
differ, the choice of methods becomes important, and this
paper has argued that the new method presented provides a
simple effective approach that should warrant consideration.
Availability
Software that implements these procedures is available electronically
on the INTERNET from the author (dtarb@
cc.usu.edu, http://www.engineering.usu.edu/dtarb/) and anonymous
ftp from fox.cee.usu.edu.
Acknowledgments. This work was supported by National Science
Foundation grant EAR-9318977. Thank you, Bill Dietrich and David
Montgomery, for use of the Tennessee Valley DEM dataset. The
DEMON code was supplied by Mariza Costa-Cabral. I would like to
thank Mariza Costa-Cabral, David Montgomery, and anonymous reviewers
for their helpful reviews.
References
Band, L. E., Topographic partition of watersheds with digital elevation
models, Water Resour. Res., 22(1), 15–24, 1986.
Beven, K. J., and M. J. Kirkby, A physically based variable contributing
area model of basin hydrology, Hydrol. Sci. Bull., 24(1), 43–69, 1979.
Beven, K., M. J. Kirkby, N. Schoﬁeld, and A. F. Tagg, Testing a
physically based ﬂood forcasting model (TOPMODEL) for three
UK catchments, J. Hydrol., 69, 119–143, 1984.
Costa-Cabral, M., and S. J. Burges, Digital elevation model networks
(DEMON): A model of ﬂow over hillslopes for computation of
contributing and dispersal areas, Water Resour. Res., 30(6), 1681–
1692, 1994.
Dietrich, W. E., C. J. Wilson, D. R. Montgomery, J. McKean, and
R. Bauer, Erosion thresholds and land surface morphology, Geology,
20, 675–679, 1992.
Dietrich, W. E., C. J. Wilson, D. R. Montgomery, and J. McKean,
Analysis of erosion thresholds, channel networks, and landscape
morphology using a digital terrain model, J. Geol., 101, 259–278, 1993.
Fairﬁeld, J., and P. Leymarie, Drainage networks from grid digital
elevation models, Water Resour. Res., 27(5), 709–717, 1991.
Freeman, T. G., Calculating catchment area with divergent ﬂow based
on a regular grid, Comput. Geosci., 17(3), 413–422, 1991.
Jenson, S. K., Applications of hydrologic information automatically
extracted from digital elevation models, Hydrol. Process., 5(1), 31–
44, 1991.
Jenson, S. K., and J. O. Domingue, Extracting topographic structure
from digital elevation data for geographic information system analysis,
Photogramm. Eng. Remote Sens., 54(11), 1593–1600, 1988.
Lammers, R. B., and L. E. Band, Automating object representation of
drainage basins, Comput. Geosci., 16(6), 787–810, 1990.
Lea, N. L., An aspect driven kinematic routing algorithm, in Overland
Flow: Hydraulics and Erosion Mechanics, edited by A. J. Parsons and
A. D. Abrahams, Chapman & Hall, New York, 1992.
Mark, D. M., Network models in geomorphology, in Modelling in
Geomorphological Systems, chap. 4, edited by M. G. Anderson, pp.
73–97, John Wiley, New York, 1988.
Marks, D., J. Dozier, and J. Frew, Automated basin delineation from
digital elevation data, Geo. Processing, 2, 299–311, 1984.
Martz, L. W., and J. Garbrecht, Numerical deﬁnition of drainage
network and subcatchment areas from digital elevation models,
Comput. Geosci., 18(6), 747–761, 1992.
Montgomery, D. R., and W. E. Dietrich, A physically based model for
the topographic control on shallow landsliding, Water Resour. Res.,
30(4), 1153–1171, 1994.
Montgomery, D. R., and E. Foufoula-Georgiou, Channel network
source representation using digital elevation models, Water Resour.
Res., 29(12), 3925–3934, 1993.
Moore, I. D., R. B. Grayson, and A. R. Ladson, Digital terrain modelling:
A review of hydrological, geomorphological, and biological
applications, Hydrol. Process., 5(1), 3–30, 1991.
Morris, D. G., and R. G. Heerdegen, Automatically drained catchment
boundaries and channel networks and their hydrological applications,
Geomophology, 1, 131–141, 1988.
O’Callaghan, J. F., and D. M. Mark, The extraction of drainage networks
from digital elevation data, Comput. Vision Graphics Image
Process., 28, 328–344, 1984.
Quinn, P., K. Beven, P. Chevallier, and O. Planchon, The prediction of
hillslope ﬂow paths for distributed hydrological modeling using digital
terrain models, Hydrol. Proc., 5, 59–80, 1991.
Quinn, P. F., K. J. Beven, and R. Lamb, The ln (a/tan ␤) index: How
to calculate it and how to use it within the Topmodel framework,
Hydrol. Process., 9, 161–182, 1995.
Tarboton, D. G., The analysis of river basins and channel networks
using digital terrain data, Sc.D. thesis, Dep. of Civ. Eng., Mass. Inst.
of Technol., Cambridge, 1989. (Also available as Tech. Rep. 326,
Ralph M. Parsons Lab. for Water Resour. and Hydrodyn., Dep. of
Civ. Eng., Mass. Inst. of Technol., Cambridge, 1989).
Tarboton, D. G., R. L. Bras, and I. Rodriguez-Iturbe, The fractal
nature of river networks, Water Resour. Res., 24(8), 1317–1322, 1988.
Tarboton, D. G., R. L. Bras, and I. Rodriguez-Iturbe, On the extraction
of channel networks from digital elevation data, Hydrol. Process.,
5(1), 81–100, 1991.
Tarboton, D. G., R. L. Bras, and I. Rodriguez-Iturbe, A physical basis
for drainage density, Geomorphology, 5(1/2), 59–76, 1992.
Wolock, D. M., and G. J. McCabe, Comparison of single and multiple
ﬂow direction algorithms for computing topographic parameters,
Water Resour. Res., 31(5), 1315–1324, 1995.
Wood, E., F. M. Sivapalan, and K. Beven, Similarity and scale in
catchment storm response, Rev. Geophys., 28(1), 1–18, 1990.
Wu, W., Distributed slope stability analysis in steep, forested basins,
Ph.D. thesis, Watershed Sci., Utah State Univ., Logan, 1993.
Zhang, W., and D. R. Montgomery, Digital elevation model grid size,
landscape representation, and hydrologic simulations, Water Resour.
Res., 30(4), 1019–1028, 1994.
D. G. Tarboton, Department of Civil and Environmental Engineering,
Utah State University, Logan, UT 84322-4110. (e-mail:
dtarb@cc.usu.edu)
(Received January 16, 1996; revised September 30, 1996;
accepted October 9, 1996.)
319TARBOTON: DETERMINING FLOW DIRECTIONS AND UPSLOPE AREAS