A triangular form-based multiple flow algorithm to estimate overland flow distribution and accumulation on a digital elevation model

perchorangeΛογισμικό & κατασκευή λογ/κού

1 Δεκ 2013 (πριν από 3 χρόνια και 6 μήνες)

72 εμφανίσεις

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.