1
A triangular form

based multiple flow algorithm to estimate overland flow
distribution and accumulation on a digital elevation model
Petter Pilesjö and Abdulghani Hasan
GIS Centre, Department of Physical Geography and Ecosystems Sciences, Lund University
S
ö
lvegatan 12, SE

223 62 Lund, Sweden
Running head: T
riangular form

based multiple flow algorithm
Corresponding Author:
Petter Pilesjö
Sölvegatan 12
223 62 LUND
Keywords
: digital terrain modelling, digital terrain analysis, hydrological modelling, surfac
e flow
esti
mation, flow routing algorithm.
2
Abstract
In t
his study
,
we present
a
newly developed
method for the
estimation
of surface flow paths on a digital
elevation model (DEM). The objective is to use a
form

based
algorithm
,
analysing flow over
single cells
by dividing them
into eight triangular facets
and to
estimate the surface flow paths on a raster DEM.
For
each cell on a
grid
ded
DEM
,
the triangular form

based multiple flow algorithm (TFM)
was used to
distribute flow to one or more of the eig
ht neighbour cells
, w
hich
determined
the flow paths over the
DEM
.
Because
each
of the eight facets covering a cell
has a constant slope and aspect, the estimations
of
—
for example
—
flow direction and divergence/convergence are
more intuitive and
less complicated
compared to
many
traditional raster

based solutions. Experiments were undertaken by estimating the
specific catchment area (SCA) over a number of mathematical surfaces, as well as on a real

world DEM.
Comparisons were made between the deri
ved SCA by the
TFM
algorithm with
eight other
algorithms
reported in the literature. The results
show that the TFM
algorithm produced the closest outcomes to the
theoretical values of the SCA compared with other
algorithms
derived
more consistent outcomes
and
was less
influenced by surface shapes.
The real

world DEM test
shows that the
TFM
was capable of
modelling flow distribution without noticeable ‘art
e
facts’, and its ability of tracking flow paths makes it
an appropriate platform for dynamic surface flo
w simulation.
3
1.
Introduction
To
pography is critical for modelling distributed hydrological processes
, and especially surface
/overland
flow
. The key param
eter in catchment topography
that has to be estimated
is
flow distribution, which
indicates how overland flow is distributed over the
terrain
. Slope
direction
controls the pathway of the
overland
flow and also substantially influences the sub

surface flow pattern.
Many authors have discussed different
flow rout
ing
algorithm
s
and how they have been applied on
digital elevation models (DEMs).
The following paragraphs are based
mainly
on
discussions by
Zhou et
al. (2
011) and Pilesjö (2008
).
One can conclude that t
he use of a
DEM
has made it possible to estimate
flow at each location over a surface. Based on the flow distribution estimation, t
he drainage pattern over
a surface
—
as well as various hydrological parameters
,
such as catchment area and up

stream flow
accumulation
—
can be
modelled
(
Wilson et al., 2008
)
.
The most critical spatial data required for digital hydrological modelling are the surface elevations,
which are typically modelled as discrete samples over the land surface. Although there have been
numerous models
developed for this purpose
—
such as the triangulated irregular network (TIN),
and
digital contour and hybrid structures (i.e.
,
a grid with intermeshed break

lines)
(
Ackerman and Krauss,
2004
)
—
the grid
ded
DEM
is
the most commonly used data source for terrain analysis
because
of its
simple data structure, ease of implementation
,
and rapidly growing applications of digital
photogrammetry and remote sensing technology
(
Li et al., 2005
)
.
Grid DEM

based hydrological m
odelling algorithms have been developed since the early stages of
geographical information systems (GIS)
(
Beven and Moore, 1994; Wilson and Gallant, 2000; Zhou et
al., 2008
)
. However
,
because
a grid DEM itself is an approximation
of
the real

world continuous surface
using regularly spaced samples, the implementation of terrain analysis models is inevitably affected by
the DEM
grid structure. In practice, numerous assumptions and optim
is
ations about mass transportation
4
and movement on
a specific local surface have to be made
to
establish a workable mathematical model
(
Holmgren, 1994
)
, resulting in significantly different approaches and methodologies towards terrain
modelling
(
Zhou and Liu, 2002
)
. It has been recogn
is
ed that such modelli
ng algorithms may produce
significantly inconsistent results related to the terrain complexity and DEM data property
(
Zhou and Liu
,
2004
; Zhou et al., 2006
)
, which lead to the presence of some significant artificial patterns (known as
‘
art
e
facts
’
) in the o
utput.
Generally speaking, on a grid DEM, the surface flow and local catchment area are approximated by
applying one out of two common
approaches, namely, the use of single flow direction (SFD)
and/
or
multiple flow direction (MFD) algorithms, according to whether
the algorithms
consider flow
divergence
s
. The SFD algorithms
(
O’Callaghan and Mark, 1984; Fairfield and Leymarie, 1991; Lea,
1992; Orlandini et al. 2003
)
do not allow flow divergence and res
trict the mass movement to
only one
downhill direction
at
a time, even though the flow may be proportionally distributed into two adjacent
downstream cells
(
e.g.
,
Tarboton, 1997
)
. The MFD algorithms
(
Freeman, 1991; Quinn et al., 1991;
Costal

Cabral and Bur
ges, 1994; Qin et al., 2007; Seibert and McGlynn, 2007
;
Gallant et al., 2011
)
consider flow divergence and assume that mass (or flow) on a grid DEM can be transported to more than
one flow direction. In many cases, SFD algorithms can produce satisfactory results over horizontally
concave surfaces where convergence of flow is assume
d. On plane or horizontally convex slopes
,
where
parallel or divergent flows are more likely, it is often more appropriate to divide the flow into two or
more directions. Combinations of the two types of algorithms are often preferred when modelling water
flow over natural surfaces
(
Pilesjö et al., 1998
)
.
Studies show that a large part of the uncertainty in derived hydrological parameters
can be explained by
the grid data structure of the DEM
(
Zhou and Liu
, 2004)
. When a natural continuous surface is
repres
ented by regularly distributed spot heights, it is inevitably difficult to determine the way that
surface flow distributes over adjacent cells
(
Olsson and Pilesjö, 2002
)
. On the other hand, flow
5
estimation over a surface with a constant slope and aspect, s
uch as a facet in a TIN, would be
considerably more consistent
,
and an SFD algorithm would be adequate without the complication of
flow divergence.
In this study
,
we
are developing
a new
flow routing
method, taking advantage of the SFD algorithms
concerning consistency, the MFD algorithm concerning the flexibility/flow diversity, the TIN modelling
regarding the continuity
, and applying it
to the grid DEM. The grid cells in the DEM
(each of them
treated as
a
‘
centre cell
’
when estimating flow from t
hat cell)
are further partitioned into triangular
facets, which have constant
slopes and aspects. Based on these
conclusions
, the surface flow path from
each facet can be uniquely tracked
and distributed to
other facets covering the centre cell,
and/
or to
one
or more of the eight cells surrounding
the
centre cell.
A comparison between the proposed algorithm and
other algorithms reported in the literature is made using
a
data

independent test method
(
Zhou and Liu
,
2002
)
, as well as
a
test on a real

world DEM
.
To summar
is
e, t
he aim of this study is to create
and evaluate
a
flow

algorithm that can simulate
overland
flow in a realistic way on a digital e
levation model, and overcome
the previous
ly
mentioned problems.
How flat areas and man

made barriers are taken care of, discussed by
Hasan
et al. (2012
a
), is
beyond the
scope of this study.
2.
Methodology
A key challenge when estimating flow over grid

based DEMs is how to model the flow movement
over
each grid cell
(
Zhou et al., 2006
)
. A typical assumption for this has been to assign a ‘flow package’ (or a
package of water) to the centre of each grid cell. Based on this, the grid cell is treated as a point on a
continuous surface
, where subsequently, a
unit of flow pa
ckage is generated and flows to the next
downhill point(s). This assumption makes it difficult to incorporate the form of the surface
overlapping
the cell
into the
flow routing
.
6
In this study, we attempt to create a computational model that is capable of s
imulating the surface flow
,
with less impact
than the
grid data structure of a DEM. The realistic simulation of the flow pattern is
considered valuable not only as a means to estimate drainage areas
but also
in estimating
the amount of
water in time and sp
ace
(
Pilesjö, 2008
)
.
The
approach is to sub

partition the grid cell according to the
surface form so that the flow direction over a local surface can be uniquely defined
(Tarboton, 1997)
.
A
grid cell is further partitioned by creating
eight local
,
triangular
facets between the cell
centre
and the
eight surrounding cell centres of the neighbouring cells
.
2.
1
Computation of slope and aspect for each triangular facet
When a facet is defined, the slope and aspect values are constant for this surface. Th
e coordinates of the
three vertices of a triangle (compare to Figure 1) are known as
C
1
(
x
1
,
y
1
,
z
1
),
C
2
(
x
2
,
y
2
,
z
2
) and
M
(
x
3
,
y
3
,
z
3
). The facet is formed as a plane
as follows
:
(1)
where
a
,
b
and
c
can be derived as:
(2)
Let
p
and
q
denote the gradients in the W

E and N

S directions, respectively. According to Equation 1,
for a triangular facet with known vertices, we have
, and
(3)
The slope (
) and aspect (
) of the facet can therefore be derived as:
c
by
ax
y
x
f
z
)
,
(
1
1
1
2
1
3
1
3
1
2
1
2
1
3
1
3
1
2
1
2
1
3
1
3
1
2
1
3
1
2
1
2
1
3
1
)
)(
(
)
)(
(
)
)(
(
)
)(
(
)
)(
(
)
)(
(
)
)(
(
)
)(
(
by
ax
z
c
y
y
x
x
y
y
x
x
z
z
x
x
z
z
x
x
b
y
y
x
x
y
y
x
x
z
z
y
y
z
z
y
y
a
a
x
f
f
p
x
b
y
f
f
q
y
7
a
a
a
b
p
p
p
q
b
a
q
p
90
arctan
180
90
arctan
180
arctan
arctan
2
2
2
2
(4)
The aspect angle is calculated
clockwise from North
.
2.
2
Flow routing within a cell
The proposed
triangular form

based multiple flow
algorithm
(TFM)
combines the advantages of
different flow distribution algorithms in a very simple but innovative way.
The work is partly based on,
and
can be seen as a further development of, the TFN algorithm presented by Zhou et al. (2011)
. The
major improvement
consider
s
form. The TFM
algorithm is
based on
multiple flow distribution allowing
overland
flow to all lower cells
surrounding a centre cell
w
hile TFN is a ‘single flow algorithm’
.
BY
dividing cells into sub

cells the proposed TFM algorithm topographic form is treated in more detail,
allowing flow to all neighbouring cells.
It is developed to be consisten
t
for all terrain types:
convex,
concave, plane
terrain, as well as t
heir combinations
.
Around the midpoint (M
, see
Figure
1
) of the cell in question
(the centre cell from where the flow is
estimated)
, eight planar triangular facets are constructed with midpoints of two adjacent c
ells (C1 and
C2). With the aid of these eight triangular facet
s, our current grid cell (centre
cell) is divided
into eight
triangular facets (
1

8 in
Figure
1
). The
slope and
slope direction
(aspect)
of each of these triangular
facets can be calculated. The
form of the current grid cell
is represented by the combined surface of the
eight triangular facets
. The area of each
facet
is
equal to 1/8
of
the cell size
,
and represents the flow
portion contributed by that facet. The
overland flow of each
triangular f
acet is to be routed
out of the
facet (
towards other facet
(
s
) or neighbouring cell(s))
,
or stays in the same triangular facet
depending on
slope direction
.
8
In the following section
s
, we
take triangular facet
number one in
Figure
1
as an example to explain the
possible
estimations of
flow routing
within
the eight facets
.
Naturally, t
he same
approach
is applied
when estimating the
flow
from another of the eight facets (2

8) to its
neighbouring
facets.
The first step is to calculate the slope
direction/
aspect value for
the facet to be
analysed
.
As
illustrated
in
Figure
2, depending on the aspect value
,
three
flow routing
alternatives are possible;
w
ater can
be
directly routed from a facet to a
neighbouring
cell. This alternative results
in
NO routing to other facets
and is
denoted
as
stay
; ALL water on a facet can be routed to ONE other
neighbouring
facet
and is
denoted
as
move
;
w
ater on a facet can be routed to a
neighbouring
facet AND a
neighbouring
cell, or to
two
neighbouring
facets
and is
denoted
as
split
.
The three alternatives are described in further detail
below:
1.
If the aspect
value
of facet 1
is between 0 and 45 degrees (see
C
ase 1,
Figure
2a
,
where
0 ≤ ASP ≤
45), the
n the
amount of water (directly proportionate to the area of the facet) to be transported from
this facet to a
neighbouring
cell will
stay
as it is (1/8 of a cell)
.
2.
If the aspect
value is between 90 and 180 degrees (see
C
ase 3,
Figure
2
d
, where
90 ≤ ASP ≤ 180)
,
or
between 225 and 270 degrees (see
C
ase 3,
Figure
2e
, where
225 ≤ ASP ≤ 270), then
all
the
water on
the facet will be
moved
to one adjacent facet
; that is, to
facet n
umber 2 and
/or
facet number 8
.
This
means that the amount of water to be routed from fac
et 1 to a
neighbouring
cell in these cases will be
0, and this water (1/8) will instead be added to facet
number
2 and
8, respectively
.
3.
If the aspect
value is between 45 and 90 degrees (see
C
ase 2,
Figure
2b, where
45 < ASP < 90)
,
or
between 270 and 360 degrees (see
C
ase 2,
Figure
2
c
, where
270 < ASP < 360), then the
water on
facet 1 will be
split
,
and partly routed into one
neighbouring
facets according to a vector split (e.g.
,
Pilesjö, 200
8
). The amount of water to
be
later routed to a
neighbouring
cell will initially stay on the
facet, while the other portion will be moved to a
neighbouring
facet. In the example
given in
Figure
9
2b
,
some water will stay (for further routing to the cell above the facet
;
see below), while some water
will be moved to facet number 2. In the example
given
in
Figure
2
c
,
some of the water will be
moved to facet number 8.
If the aspect
value is between 180 and 225 degrees (see
C
ase 4,
Figure
2
f
,
where 1
80 < ASP < 225), then the
water on facet 1 will
be
split
between two adjacent facets;
one
portion will be added to
the
amount of water (proportional to the area)
of facet n
umber
2
,
while
another
portion will be added to facet n
umber
8. The
resulting amount of water of facet number
1
will be 0.
Note, in the case of two facets sloping towards each other, i.e., contradictive slopes, no water is routed
between the two facets, but the alternative
stay
is applied.
Because water can be transported between more than two facets (from facet 1 to facet 2; to facet 3 to
facet 4; and then to neighbouring cells), flow routing between facets continues until all water has
reached the
outflow
facet(s)
of the cell
. Water is the
n transported to neighbouring cells. The result of
this first step is the redistribution of water within a cell. It also estimates where (to which neighbour
cells) water will route.
2.3 Flow distribution
to
neighbouring
cells
After doing the
flow routing
for
each cell
—
which
consist
s
of eight individual facets
—
the routing of
water to adjacent cell(s) takes place. In one or more of the eight facets
,
water is
‘
waiting
’
to be routed to
neighbouring
cell(s). If the cell represents a
concave
landform
,
we can hav
e only one facet holding all
the water (directly proportional to the area, 1, of one cell)
;
if the cell represents a convex surface
(compare
to
a
pyramid)
,
there might be water in all eight facets.
Each facet has two
neighbouring
cells (e.g.
,
see
cells C1
and C2 for facet number one in
Figure
1)
,
and
the next step is to distribute the
accumulated water
between these cells.
Two different cases can
then
be
found
:
10
1.
If one
neighbour
ing
cell is lower
in elevation
than the
centre
cell
(where the facet is located)
and the
other cell is higher or equal
to the
centre
cell,
then the
water
accumulated in the
facet will all be
distributed to the
lower cell.
2.
If both
neighbouring
cells are lower in elevation than the
centre
cell, then
the
accumulated
water in
the facet
will be
proportionally
distributed to both
lower
cells
to slope (
an
x
value of 1
,
Equation
5
)
:
where
i,j
= flow directions (1…
2
) to lower cells,
f
i
= flow proportion (0…1) in direction
i
, tan
ß
i
=
slope gradient between the centre cell and the cell in direction
i
, and
x
= variable exponent.
3.
If both neighboring cells are higher
or equal
in elevation than the central cell and there are other
lower elevation cells then the accumulated water in the facet will be proportionally distributed to
all lower cells.
2.4 Accumulated flow estimation and specific catchment area
Flow distribution to
ne
ighbouring
cells is estimated for all cells except for the border cells in the DEM.
Starting f
rom
cells with no incoming water
—
that is,
peaks or cells at the edge of the DEM
—
flow
accumulation
s
(in the unit of cells)
are
estimated
as going
downhill in the catchments.
The result
s are
estimated
as
values of flow accumulation for every cell
,
with the exception of
the border cells in the
DEM.
When comparing results between different flow routing algorithms
,
the specific catchment area (SCA)
was
used.
The SCA can be estimated using
the definition and method reported by Costa

Cabral and
Burges
(
1994
)
as:
2
1
tan
tan
j
x
j
x
i
i
f
1.
, for all ß > 0
(5)
11
)
cos
sin
(
A
A
g
TCA
TCA
SCA
(6)
where
SCA
is the specific catchment area defined as the upstream catchment area per unit contour line
(m
2
/m);
TCA
is the total upstream catchment area at a given point (m
2
);
is the flow width of the
catchment outlet (m);
A
is the slope aspect at the centre of the
grid cell (º)
, estimated by approximation
of a second order trend surface and
representing the general flow direction of that cell; and
g
is the grid
size (m).
3. Experiment
Two different methods have been used to test the proposed flow estimation algorithm (TFM).
The first
one is a
data

independent test based on the method proposed by Zhou and Liu
(
2002
)
, testing the results
of flow estimation on a number of mathematical surf
aces. This makes it possible to derive and compare
quantitative measurements o
f
the accuracy for different methods.
We also apply the proposed algorithm
to a real

world DEM. The spatial pattern
s
of estimated flow and accumulated flow
are
then visually
exam
ined
,
to
detect if there are significant
artefacts
.
3.1
Data

independent assessment
Repeating the experiment carried out by Zhou et al. (2011)
,
the following surfaces were used: ellipsoid
(representing a convex slope), inverse ellipsoid (representing a con
cave slope), plane (representing a
straight slope) and saddle (representing
a
convex/concave slope association). The four different
mathematical surfaces are visualised
below,
in
Figure
3
.
The generation of the surfaces, including
the
formula
defining the
elevation at each and every cell
,
can be found in Zhou et al. (2011).
According to
Zhou et al. (2011), the
cell size for all surfaces
is
5
metre
s
, and the number of ce
lls varies between 15
000 and 40
000. For the plane surface, the slope is set to
66º,
and
the
aspect
to
243º.
12
At each cell of the mathematical surfaces
,
the
theoretical ‘true’ specific catchment area (SCA)
was
calculated according to the method
described by Zhou and Liu
(
2002
)
.
When the ‘true’ SCA values
,
as
well as the estimated flow algorithm dependant
ones are known, the
error at each cell can be computed
as:
t
i
i
SCA
SCA
E
(7
)
where
E
i
denotes the error (or residual) at the
i
th cell using a selected algorithm,
SCA
t
and
SCA
i
denote
the theoretical and
estimated
value (using the
flow
algorithm
to be tested
) of the SCA at the
i
th cell,
respectively. The Root

Mean

Standard Error (RMSE), Mean Error (ME, denoted as
E
) and Standard
Deviation of the residuals (SD, denoted as
) were then computed
according to Zhou et al. (2011)
for
the assessment and comparison of the selected algorithms:
n
E
RMSE
n
i
i
1
2
(8
)
n
E
E
n
i
i
1
(9
)
1
1
2
n
E
E
n
i
i
(10
)
where
n
denotes the total number of grid cells for each mathematical surface DEM.
The proposed
TFM
algorithm was applied
on the four mathematical surfaces
,
and the estimated flow
accumulation values for each cell were converted to SCA in the unit of grid cells ac
cording to
Equation
6.
Then,
the estimated SCA values were compared with the ‘true’ values computed. Statistical analyses
13
of the differences (RMSE,
E
, and
) were carried out
to
estimate overall accuracy as well as spatial
distribution of errors.
3.2 Selection of algorithms for comparison
The proposed TFM algorithm was compared with eight commonly used, tested, and
well documented
flow algorithms. T
hree
of these can be classified
as grid based single flow algorithms (D8,
D8

LTD
,
and
D
),
three
as grid based multiple flow algorithms
(MD
, FMFD, and QMFD)
, and two as (partly)
vector

based algorithms
(DEMON and TFN)
. The selected algorithms are
as follows
:
Deterministic Eight

Node
(D8)
(
O’Callaghan and Mark
, 198
4
)
Deterministic Eight

Node Least Transversal Deviation
(D8

LTD)
(
Orlandini et al., 2003; Orlandini
and Moretti, 2009
)
Deterministic
infinite number of possible single flow directions
(D
)
(
Tarboton, 1997
)
Triangular
Multiple Flow Direction
(MD
)
(
Seibert and McGlynn, 2007
)
Freeman Multiple Flow Direction
(FMFD)
(
Freeman, 1991
)
Quinn Multiple Flow Direction
(QMFD)
(
Quinn et al., 1991
)
Digital Elevation Model Networks
(DEMON)
(
Costa

Cabral and Burges, 1994
)
T
riangular Facet Network
(TFN) (Zhou et al., 2011)
W
hich algorithm to select for comparisons
can
always be discussed
. However, the selected eight
algorithms
were judge
d
to represent a wide variety of solutions, involving raster as well as vector

based
algor
ithms, all
widely available to the research community. W
e used
the
System for Automated Geo

Scientific Analysis
,
SAGA
(2012)
for
testing
D8, D
, MD
, FMFD and DEMON algorithms, the
14
original FORTRAN code available from the author’s website
1
for the D8

LTD algorithm
, and the
original code for the TFN algorithm
available
directly from the authors
. Software implementation for the
QMFD algorithm was conducted using the C++ programming language
from the Zhou et al. study
(2011)
. ArcGIS
2
GIS softw
are was used for comparison and visual
is
ation of the results.
3.3 Real

world testing
To
visually examine the result of the proposed TFM algorithm
,
it was applied
to
a real

world DEM. This
1
metre resolution DEM was interpolated form LiDAR data measured at
the Stordalen mire and its
catchment area. Stordalen is a peat
land area in the Arctic region 10 km west of Abisko (68º 20' N, 19º
03' E), north of Sweden. The LiDAR data
,
as well as the interpolation procedure
,
are described in detail
in
Hasan
et al.
(2012
b
).
Because
the differences between the different flow algorithms are relatively small (see below), the main
purpose of real

world t
esting is to identify possible
artefacts
created by the proposed TFM algorithm.
However,
to
highlight
differences,
we h
ave also
included
accumulated flow estimations based on the D8
algorithm.
4. Results
The tests of the different flow algorithm can be divided into two parts:
first, the
quantitative test of the
significant difference
s between estimated values and ‘true values’ of SCA over the mathematical
surfaces, and
second, the
qualitative test examining spatial distribution of errors over the mathematical
surfaces, as
well as visual interpretation of the real

world
modelling
.
4.1
Quantitative accuracy assessment
1
http://www.idrologia.unimore.it/orlandini/download.html
.
2
ArcGIS 9.3, © Environmental Systems Research Institute, Inc., http://www.esri.com
.
15
The
significant difference
s between the
nine
different flow algorithms and the ‘true values’ of SCA for
the four mathematical surfaces are presented in Table 1.
For three out of four of the mathematical
surfaces, t
he result
s show that the proposed TFM algorithm
outperforms the other
eight
algorithms in
terms of RMSE. It has the lowest RMSE values
out of
5.88
m
, 8.28
m
, and 7.66
m
, for the ellipsoid, the
inverse ellipsoid
,
and the saddle
, respectively.
For the plane, the
TFM
algorithm has a higher (6.13
m
)
RMSE value than the TFN algorithm (1.86
m
)
and is
also just higher than the DEMON algorithm (5.18
m
).
Regarding the systematic bias, represented by the Mean Error (
E
) value
,
the TFM algorithm shows
better results than the
all
other algorithm for two of the four mathematical surfaces: the ellipsoid (
E
=

0.99
m
) and the saddle (
E
= 0.61
m
). For the inverse ellipsoid the differences
between the TFM result (
E
=
4.00
m
) and the slightly better TFN (
E
=

3.95
m
) and DEMON (
E
=
3.88
m
) are marginal.
For
the plane, only the TFN algorithm shows a better result (
E
= 1.85
m
) than the TFM algorithm (
E
=
3.12
m
).
If we assume normal distribution, the
n the
Standard Deviation of the Residuals (
) indicates
the
representability of the
M
ean
E
rror; the lower
the S
tandard
D
eviation, the more narrow distribution of
errors around the
M
ean
E
rror.
For the ellipsoid surface, the TFM algorithm has a
S
tandard
D
eviation of
5.80 m, which is the third lowest of the eight algorithms. However
, because
the FMFD (
= 2.32 m) and
the QMFD
(
= 4.61 m) algorithm both have considerably higher mean errors (6.87 m and 6.64 m
compared to

0.99 m) the deviation from zero (no error) is actually smaller for the TFM algorithm. This
is also the case for the saddle, where only the FMFD algorithm has a
lower
S
tandard
D
eviation than
TFM (6.95 m compared to 7.63 m)
,
but a much higher mean error (13.19 m compared to 0.61 m).
For
the plane
,
the TF
M
standard deviation (5.27 nm) is higher than the
deviations of
DEMON (3.91 m) and
16
TFN (0.11 m) algorithms, and
f
or the inverse ellipsoid
the TF
M
algorithm shows the lowest standard
deviation
out
of all
the
tested algorithm
s
.
The
M
ean
E
rror and
S
tan
dard
D
eviation also indicate possible
skewness
in the estimation of flow
accumulation. If the standard deviation is
less than the absolute value of the mean error
,
this
indicates
that at least 83% of the errors are either all positive or all negative. Such a systematic error is not
desirable.
Referring to
Table
1, the absolute mean error never exceeded the standard devi
ation for t
hree
algorithms (DEMON
, TFN,
and TFM)
;
exceeded it
at
one surface
,
for
three
algorithms (D8, D8

LTD,
and
D
)
;
exceeded it
at
two surfaces
,
for one algorithm (MD
)
;
and exceeded it
at
three out of four
surfaces
,
for two algorithms (FMFD and
QMFD).
When comparing the frequency distribution of
positive and negative errors for these eight algorithms, only the TFM algorithm did not exceed skewness
of the proportion 80

20 (%) for any of the four surfaces.
D8

LTD exceeded 80

20 for one surface, D8,
D
, MD
, and TFN did it for two surfaces, DEMON for three, and FMFD and QMFD exceeded the 80

20 proportion for all four surfaces. It is also worth noting that the skewness
was
one

sided (either over
or underestimating SCA for all mathematical surfaces) fo
r all tested algorithm but TF
M
.
In
Figure
4
we illustrate the results of one single flow algorithm (D8), one multiple flow algorithm
(
QMFD
), and the proposed form

based algorithm (TFM) applied on the four mathematical surfaces.
Differences in flow routing
are obvious, especially for the saddle
, where only the T
F
M algorithm
logically routes 50% of the water in NW and SE directions
.
Regarding the quantitative comparison between the proposed triangular form

based multiple flow
algorithm (TFM) and the eight alt
ernative algorithms it clearly shows that the TFM algorithm produces
the best over

all results. Even if it did not produce the best result in every test over every mathematical
surface, it shows by far the best and most consistent results. Referring to Tab
le 1 one can conclude that
out of the 12 measurements of accuracy (RMSE,
E
, and
for the four mathematical surfaces) the TFM
17
algorithm outperforms all other algorithms in six of these. The algorithms producing the second best
result,
the TFN algorithm, outperforms all other algorithms in three out of these 12 cases. It should also
be noted that these three cases are all linked to estimation on the plane surface.
Generally it can be concluded that the raster

based single flow algorithms
(D8 and D8

LTD) produce
the poorest results, with RMSE values between 50 and 224 m (average = 130.0 m). The raster

based
multiple flow algorithms (D
, MD
, FMFD, and QMFD) produce better results, with RMSE between 7
and 66 m (average = 32.7 m). The (partl
y) vector

based algorithms (DEMON and TFN) produce even
more reliable results, with RMSE between 2 and 56 (average = 14.7 m). The reason for this can be
explained by the obvious over

simplification adopted by the single flow algorithms, in combination with
the advantages in geometric precision connected to the vector solutions. The average RMSE value for
the TFM algorithm is 10.0 m.
4.2 Spatial distribution of errors on the artificial surfaces
The ‘true values’ of SCA
,
as well as the estimated values applying the proposed TFM algorithm for the
four mathematical surfaces
,
are presented in
Figure
5. Even if differences can be
observed
, it is clear that
the modelled SCA values
,
to a high degree
,
coincide with the true
valu
es
.
T
he spatial distribution of residuals of derived SCA
—
with respect to the theoretical ‘true’ value of SCA
on the four mathematical surfaces
,
for all
of the
nine tested algorithms
—
is shown in
Figure
6
.
Logically
,
and
based on the errors reported above (see
Table
1
, for example
)
, the TFM, TFN and DEMON
algorithms
show relatively less errors compared to the other algorithms tested. A general difference
between the raster

based algorithms (
D8, D8

LTD,
D
, M
D
, FMFD, and
QMFD) and the vector

based
algorithms
(DEMON and TFN)
can also be identified
:
t
he vector

based
algorithms show a more
random pattern of error distribution,
while the errors in the raster

based
algorithms seem to be more
additive. However, for the TFM algor
ithm, this systematic bias seems to be relatively low.
18
The less

biased pattern
of the TFN algorithm is explained by the way SCA values are calculated. By
counting the number of flow lines passing each cell
,
in vector mode
,
the derived SCA values of
the
upper stream cells do not affect that of the down

stream cell
(see Zhou et al., 2011)
.
However, as will be
further discussed, the TFN algorithm do
es
not produce better
overall
results
and is
much more CPU
demanding.
4.3 Real

world testing
Even if visual
differences are often difficult to detect when modelling flow over real

world surfaces
,
the
derived SCA values between the
TFM
and D8 algorithm w
ere
estimated for a peat land area in northern
Sweden.
Figure
7
presents the visual comparison between the
resu
lts of the two algorithms.
The results
show that the
TFM
has produced a more balanced simulation on the flow distribution on this terrain
surface, with a more realistic spatial pattern on the
different terrain types (
straight slope
s, converging
slopes and
diverging slopes).
In comparison, the D8 algorithm show
s a possible advantage in
representing the drainage network, but it also demonstrates its mostly critic
is
ed weakness
:
an
unacceptable 45

degree bias that results in significant
artefacts
,
such as
unrealistic
parallel drainage
li
nes
.
5. Discussion
Even if the proposed TFM algorithm takes a raster approach when estimating the overland flow
pattern
—
although this
has been recogn
is
ed as one of the major sources of error for flow estimation
(
Zhou and Liu
, 2004)
—
it seems superior to all other tested algorithms, including the (partly
)
vector
based DEMON and TFN algorithms.
5.1 Accuracy comparison
19
The accuracy estimations have been carried out using
both
quantitative comparisons between estimated
flow accumu
lation/SCA and ‘true values’ over mathematical surfaces, and visual interpretation of a
natural surface in northern Sweden. The main purpose of th
e latter is to reveal possible
artefacts
,
caused
by over

simplifications and ‘illogical flow routing’.
C
ommon
problem
s
connected to the multiple flow algorithms
include
the limited possibilities in
adjusting the overland flow to
the
terrain form. Logically
,
we assume the flow to convert over concave
surfaces and divert over convex surfaces. This can be regulated b
y the
x
exponent in
Equation
5, but this
is rarel
y
performed
(see Pilesjö et al., 1998
). However, this is taken care of in the proposed form

based
algorithm. The logic in the flow routing is exemplified in
Figure
4, where it can be
observed
how water
is an
d should be divided over a ridge (the saddle).
T
he results of visual analysis over both artificial and real

world surfaces show no obvious systematic
bias or artificial spatial pattern
for the TFM algorithm
.
The single flow algorithms have the capability of
modelling concave landforms relatively well, but show illogical results for plane and convex
landforms
.
The opposite is valid for most of the multiple flow algorithms;
because
flow is normally distributed t
o
all neighbour cells with lower elevation this creates unrealistic flow accumulation values for concave
surfaces.
5.2 Potential implementation
The source code of the proposed TFM algorithm is available from the authors. In this code, also
solutions how to
tackle
e.g.,
sinks, flat areas, and man

made barriers
are implemented
(Hasan et al
.
,
2012a).
Compa
red to
vector based methods (
e.g.,
the TFN algorithm), the computational time (CPU) of
the proposed
TFM
algorithm is also fast
, making it more useful
. Most
likely
, the major hurdle for the
implementation of the TFM algorithm, as well as other more advanced multiple flow algorithms, is the
slow response from the research and industry communities.
This, in combination with the question if
20
‘better’ flow algorith
ms are really needed (and they not always are), will probable slow down the
potential implementation.
5.3 Assessment methods
The proposed TFM method has been evaluated and compared with eight other commonly used
algorithms reported in the literature. Both
quantitative and qualitative evaluations of flow
accumulation/SCA, on mathematical surfaces (ellipsoid, inverse ellipsoid, plane, and saddle) as well as
real

world data, have been carried out. Root mean standard errors, mean errors, and standard deviations
of the residuals have been calculated and compared. Visual interpretations have been made.
However,
m
ore thorough and rigorous evaluation methods
can always be
performed
,
not only on the quantitative
error
s
of the estimated catchment area (
such as
the
area estimation error)
but also
on the uncertainty of
the flow path or their derivatives’ morphology (
see
Orlandini and Moretti, 2009
)
.
6. Conclusion
In this study we have tested a newly developed triangular form

based multiple flow algorithm (TFM).
The al
gorithm is form

based in the
sense
that it divides
the cell from where the flow routing takes place
into eight triangular facets
to
estimate flow
over
the cell. The flow is
then
distributed to
one or more of
the eight
neighbour
cells, determining the flow paths over the DEM.
Because
each of the eight facets
covering a cell has a constant slope and aspect, the estimations of, for example, flow direction and
divergence/convergence are
more intuitive and
less complicated compared t
o
many
traditional raster

based solutions. Experiments undertaken by estimating the specific catchment area (SCA) over a
number of mathematical surfaces, as well as on a real

world DEM
,
clearly indicate
an
all

round
advantage of estimating flow and flow ac
cumulation in all forms of slopes
when using the TFM
algorithm compared to other algorithms
. The TFM algorithm outperforms eight commonly used
compared algorithms
(D8, D8

LTD,
D
, MD
, FMFD, QMFD,
DEMON, and TFN)
, proven by both
21
quantitative measurements
including RMSE, ME and SD, and qualitative evaluation including the visual
analysis of spatial patterns of residuals on the artificial surfaces and SCA over a real

world DEM.
It is also concluded that the fact that the TFM algorithm is raster

based in term
s
o
f flow distribution
between different cells in the DEM, making it considerably faster than vector based alternatives (
e.g.,
TFN), is highly desirable when processing larger DEM.
Overall,
the proposed TFM algorithm is considered both
adequate in represen
ting overland flow over a
DEM, resulting in
accurate estimations of flow accumulation and SCA, and efficient enough to have a
high potential for broad implementation.
7.
Acknowledgements
The European Union funding programme Erasmus Mundus
‘
External
Cooperation Window
’
(EMECW
lot8) financed the second author of this study. The authors are grateful to Qiming Zhou, Hong Kong
Ba
p
tist University, and Yumin C
hen,
Wuhan University
, for fruitful cooperation and data exchange.
8. References
Ackermann, F., K.
Kraus, (2004), Reader commentary: grid based digital terrain models,
Geoinformatics, 7(6), 28

31.
Beven, K. J., and I. D. Moore, (eds.) (1994), Terrain Analysis and Distributed Modelling in Hydrology,
John Wiley & Sons, Chichester.
Costa

Cabral, M. C., and
S. J. Burges (1994), Digital elevation model networks (DEMON): A model of
flow over hillslopes for computation of contributing and dispersal areas, Water Resources Research,
30(6), 1681

1692.
22
Fairfield, J., and P. Leymarie (1991), Drainage network from gr
id digital elevation models, Water
Resources Research, 27(5), 709

717.
Freeman, T. G. (1991), Calculating catchment area with divergent flow based on a regular grid,
Computers & Geosciences, 17(3), 413

422.
Gallant, J.C., and Hutchinson, M.F., 2011. A
differential equation for specific catchment area, Water Re
sources
Research, 47(5), W05535
,
http://dx.doi.org/10.1029/2009WR008540
Hasan, A., P. Pilesjö and A. Persson (2012
a
) Drainage Area Estimation in Practice, how to tackle
artefacts in real world data, In the proceeding Surface models for geosciences, Symposi
um GIS Ostrava
2012 Ostrava, Czech republic
.
Hasan, A., Pilesjö P. and Persson A. (2012b) On generating digital elevation models from LiDAR data
–
resolution versus accuracy and topographic wetness index indices in northern peatlands, Geodesy and
Cartography, Taylor & Francis.
Holmgren, P. (1994), Multiple flow direction algorithms for runoff modelling in grid based elevation
models: An empirical evaluation, Hydrological Processes, 8, 327

334.
Lea N. J. (1992), An aspect

driven kinematic routing al
gorithm, in Overland Flow: Hydraulics and
Erosion Mechanics, edited by A. J. Parsons, and A. D. Abrahams, pp. 393

407, UCL Press, London.
Li, Z., Q. Zhu, and C. Gold (2005), Digital Terrain Modelling: Principles and Methodology, CRC Press,
Boca Raton.
O’C
allaghan, J. F. and D. M. Mark (1984), The extraction of drainage networks from digital elevation
data, Computer Vision, Graphics, and Image Processing, 28, 323

344.
23
Olsson, L. and P. Pilesjö (2002), Approaches to spatially distributed hydrological modelli
ng in a GIS
environment, in Environmental Modelling with GIS and Remote Sensing, edited by A. Skidmore,
Taylor & Frances, London.
Orlandini, S., G. Moretti, M. Franchini, B. Aldighieri and B. Testa (2003), Path

based methods for the
determination of non di
spersive drainage directions in grid

based digital elevation models, Water
Resources Research, 39(6), 1144, doi:10.1029/2002WR001639.
Orlandini, S., and G. Moretti (2009), Determination of surface flow paths from gridded elevation data,
Water Resources Res
earch, 45, W03417, doi: 10.1029/2008WR007099.
Pilesjö, P. (2008), An integrated raster

TIN surface flow algorithm, in Advances in Digital Terrain
Analysis, edited by Q. Zhou, B. Lees, and G. Tang, pp. 237

255, Springer, Berlin.
Pilesjö, P., Q. Zhou, and L.
Harrie (1998), Estimating flow distribution over Digital Elevation Models
using a Form

Based Algorithm, Geographic Information Science, 4, 44

51.
Qin, C., A.

X. Zhu, T. Pei, B., Li, C. Zhou, and L. Yang (2007), An adaptive approach to selecting a
flow

par
tition exponent for a multiple

flow

direction algorithm, International Journal of Geographical
Information Science, 21(4), 443

458.
Quinn, P., K. Beven, P. Chevallier, and O. Planchon (1991), The prediction of hillslope flow paths for
distributed hydrologi
cal modelling using digital terrain models, Hydrological Processes, 5, 9

79.
SAGA 2.0.5, © Saga User Group Association,
http://www.saga

gis.org
. Accessed 2012

04

02.
Seibert, J. and B. L McGlynn (2007), A new triangu
lar multiple flow direction algorithm for computing
upslope areas from gridded digital elevation models, Water Resources Research, 43, W04501.
24
Tarboton, D. G. (1997), A new method for the determination of flow directions and upslope areas in grid
digital
elevation models, Water Resources Research, 32(2), 309

319.
Wilson, J. P., and J. C. Gallant (eds.) (2000), Terrain Analysis: Principles and Applications, John Wiley
& Sons, New York.
Wilson, J. P., G. Aggett, Y. Deng, and C. S. Lam (2008), Water in the l
andscape: A review of
contemporary flow routing algorithms, in Advances in Digital Terrain Analysis, Edited by Q. Zhou, B.
Lees, and G. Tang, pp. 213

236, Springer, Berlin.
Zhou, Q., and X. Liu (2002), Error assessment of grid

based flow routing algorithms
used in
hydrological models, International Journal of Geographical Information Science, 16(8), 819

842.
Zhou, Q., and X. Liu (2004), Analysis of errors of derived slope and aspect related to DEM data
properties, Computers and Geosciences, 30: 369

378.
Zhou, Q., B. Lees, and G. Tang (eds.) (2008), Advances in Digital Terrain Analysis, Springer, Berlin,
462p.
Zhou, Q., X. Liu, and Y. Sun (2006), Terrain complexity and uncertainties in grid

based digital terrain
analysis, International Journal of Geographi
cal Information Science, 20(10), 1137

1147.
Zhou, Q., P. Pilesjö, and Y. Chen (2011), Estimating surface flow paths on a digital elevation model
using a triangular fa
cet network.
Water Resources Research 47 (7), art. no. W07522.
http://dx.doi.org/10.1029/2010WR009961
.
25
Table 1
. Accuracy
comparison between the TFM
algorithm and eight other selected algorithms on the
four mathematical surfaces (unit = metres). Partly from Zhou et al. (2011).
Ellipsoid
Inverse ellipsoid
Plane
Saddle
D8
RMSE
63.53
224.31
154.53
103.68
E

12.07

86.22

109.87

2.06
62.37
207.08
108.67
103.66
D8

LTD
RMSE
97.21
174.98
50.08
172.01
E

14.04

26.21

42.54
11.52
96.19
173.01
26.43
171.62
D
RMSE
21.51
66.03
62.75
27.67
E

10.65

21.43

49.92

5.71
18.69
62.45
38.02
27.07
MD
RMSE
14.52
56.37
62.75
17.32
E

11.81

18.55

49.92

4.41
8.46
53.23
38.02
16.74
FMFD
RMSE
7.25
34.94
29.62
14.91
E
6.87
14.11
25.26
13.19
2.32
31.97
15.47
6.95
QMFD
RMSE
8.08
41.24
42.26
15.76
E
6.64
17.09
31.61
13.11
4.61
37.53
28.06
8.75
DEMON
RMSE
17.21
11.21
5.18
55.64
E
4.47
3.88
3.41
6.67
16.62
10.52
3.91
55.23
TFN
RMSE
8.42
10.25
1.86
8.11
E
1.97

3.95
1.85

0.92
8.18
9.46
0.11
8.06
TFM
RMSE
5.88
8.28
6.13
7.66
E

0.99
4.00
3.12
0.61
5.80
8.23
5.27
7.63
26
Figure 1.
In a 3 by 3 cells window, the centre cell is divided into eight triangular facets (1

8). Each facet
is formed from three points; one is the centre cell (M), and the two other
are two adjacent cells (e.g.,
C1
and C2).
1
2
3
4
5
6
7
8
C1
C2
M
27
Figure 2.
An illustration of how water can be routed from one
triangular facet to other facets. Different
aspect values lead to different actions (stay, move and split
).
2.
a
1
2
3
4
2
3
3.
b
4.
c
5.
d
6.
e
7.
f
28
Figure 3
. Illustrations of the four different mathematical surfaces used in the
data

independent
assessment of flow algorithms.
Saddle
Plane
Ellipsoid
Inverse Ellipsoid
29
Figure 4
. Distribution of flow from a centre cell to
the eight neighbouring cells using different flow
algorithms over sub

surfaces (3 x 3 cells) of four mathematical surfaces (inverse ellipsoid, plane,
ellipsoid, and saddle).
Ellipsoid
Inverse ellipsoid
Plane
Saddle
100
100
99.8
98.8
99.8
100.4
99.6
100
100.7
98
100
102
101
99
97
99
101
103
93.1
100
107.2
100.4
93.1
85.9
100.4
107.2
114.2
99.3
100
99.3
98.6
100
98.6
97.3
98
97.3
100
13
55
13
19
12
12
21
26
21
4
4
34
17
37
12
29
29
42
100
100
100
7.3
35.4
7.3
50
21.1
21.1
18.1
12.8
18.1
4.4
4.4
46.1
10.6
38.4
4.9
31
31
38
Elevation
D8
QMFD
TFM
30
Figure 5
. A comparison
between the ‘true’ SCA values (A1, B1, C1, D1) and the ones modelled by
applying the proposed TFM algorithm (A2, B2, C2, D2) for four different mathematical surfaces
(ellipsoid, inverse ellipsoid, plane and saddle).
Low
SCA values
(dark brown) to high
SCA
values
(dark
blue).
A1
A2
B1
B2
D2
D1
C2
C1
31
Ellipsoid
Inverse Ellipsoid
Plane
Saddle
D8
D8

LTD
D
MD
FMFD
QMFD
DEMON
<

1,000

1,000


100

100


10

10


1

1

1
1

10
10

100
100

100
> 1,000
<

1,000

1,000


100

100


10

10


1

1

1
1

10
10

100
100

100
> 1,000
DEMON
32
TFN
TFM
Figure 6.
The spatial distribution of residuals of the derived SCA compared to the theoretical ‘true’
values over four mathematical surfaces.
33
Figure 7.
A comparison between the D8 and TFM algorithms on estimated SCA over a real

world
DEM.
TFM
D8
34
Figure captions
Figure 1.
In a 3 by 3 cells window, the centre cell is divided into eight triangular facets (1

8). Each facet
is formed from three points; one is the centre cell (M), and the two other
are two adjacent cells (e.g.,
C1
and C2).
Figure 2.
An illustration of how water
can be routed from one triangular facet to other facets. Different
aspect values lead to different actions (stay, move and split
).
Figure 3
. Illustrations of the four different mathematical surfaces used in the data

independent
assessment of flow algorith
ms.
Figure 4
. Distribution of flow from a centre cell to the eight neighbouring cells using different flow
algorithms over sub

surfaces (3 x 3 cells) of four mathematical surfaces (inverse ellipsoid, plane,
ellipsoid, and saddle).
Figure 5
. A comparison
between the ‘true’ SCA values (A1, B1, C1, D1) and the ones modelled by
applying the proposed TFM algorithm (A2, B2, C2, D2) for four different mathematical surfaces
(ellipsoid, inverse ellipsoid, plane and saddle).
Figure 6.
The spatial distribution of re
siduals of the derived SCA compared to the theoretical ‘true’
values over four mathematical surfaces.
Figure 7.
A comparison between the D8 and TFM algorithms on estimated SCA over a real

world
DEM.
Σχόλια 0
Συνδεθείτε για να κοινοποιήσετε σχόλιο