scieee Science in your language
[en] (orig)
This version is available at http://dx.doi.org/10.14279/depositonce-7231
This work is licensed under a CC BY-NC-ND 4.0 License (Creative
Commons Attribution-NonCommercial-NoDerivatives 4.0 International). For
more information see https://creativecommons.org/licenses/by-nc-nd/4.0/.
Terms of Use
Ö
zgen, I., Liang, D., & Hinkelmann, R. (2016). Shallow water equations with depth-dependent anisotropic
porosity for subgrid-scale topography. Applied Mathematical Modelling, 40(17–18), 7447–7473.
https://doi.org/10.1016/j.apm.2015.12.012
Ö
zgen, I.; Liang, D.; Hinkelmann, R.
Shallow water equations with depth-
dependent anisotropic porosity for
sub
g
rid-scale topo
g
raphy
Accepted manuscript (Postprint)Journal article |
Shallow water equations with depth-dependent anisotropic
porosity for subgrid-scale topography
Ilhan Özgen
a, ∗, Dongfang Liang
b,
Reinhard Hinkelmann
a
a
Chair of Water Resources and Modeling of Hydrosystems, Technische Universität Berlin, Germany
b
Department of Engineering, University of Cambridge, UK
a b s t r a c t
This paper derives a novel formulation of the depth-averaged shallow water equations with
anisotropic porosity for computational efficiency reasons. The aim is to run simulations on
coarser grids while maintaining an acceptable accuracy through the introduction of poros-
ity terms, which account for subgrid-scale effects. The porosity is divided into volumetric
and areal porosities, which are assigned to the cell volume and the cell edges, respectively.
The former represents the volume in the cell available to flow and the latter represents the
area available to flow over an edge, hence introducing anisotropy. The porosity terms are
variable in time in dependence of the water elevation in the cell and the cumulative dis-
tribution function of the unresolved bottom elevation. The main novelty of the equations
is the formulation of the porosities which enables full inundation of the cell. The applica-
bility of the equations is verified in five computational examples, dealing with dam break
and rainfall-runoffsimulations. Overall, good agreement between the model results and a
high-resolution reference simulation has been achieved. The computational time decreased
significantly: on average three orders of magnitude.
1. Introduction
Shallow water models are applied in a broad range of areas such as river hydraulics [1–4] , dam break simulations [5,6] ,
urban flooding [7–10] and recently also for overland flow in natural catchments [11–16] and urban runoff [17,18] , among
many others (cf. e.g. [19] ). In overland flow simulations, usually there is a large difference between the scales of the fea-
tures significantly influencing the flow and the scale of the simulation domain. For example, in a natural catchment with
a scale around a square kilometer, local depressions and microtopograpy with horizontal scales smaller than a square me-
ter influence the flow field significantly [20–22] . Similar observations are made for urban flood models where the scale of
buildings is exceeded by the scale of the city in several orders of magnitude, e.g. a building has a scale of around 100 m
2
while the city may span up to 100 km
2
. Recent developments in survey technology such as light detection and ranging (LI-
DAR) and laser scanning are able to provide high accuracy high-resolution elevation data sets at relatively low cost [23] .
However, the integration of these data into numerical models is often challenging because of finite computer resources
[24,25] . In order to capture the impact of the smallest relevant scale on the flow, the microtopography has to be explicitly
discretized. This leads to meshes with small cell size and therefore high cell number which in return leads to an increased
∗Corresponding author. Tel.: +49 3031472428.
E-mail address:
ilhan.oezgen@wahyd.tu-berlin.de (I. Özgen).
1
computational effort. Despite developments in CPU power, high-resolution simulations across large catchments are in prac-
tice often unfeasible without supercomputers [26] .
Instead of explicitly discretizing the small-scale topography, its influence on the flow can be conceptually accounted for
on coarser meshes to reduce the computational effort [27] . One such approach introduces a porosity term into the shal-
low water equations, which refers to the fraction of a computational cell available for flow and is a concept borrowed
from groundwater flow modeling. The porosity then conceptually accounts for subgrid-scale topography. In literature, the
extended shallow water equations incorporating this porosity are called shallow water equations with porosity or porous
shallow water equations. The initial porous shallow water equations have been derived by Defina [11] to account for micro-
topography in overland flow. Later, the concept has been applied in urban flood modeling as a building treatment method
[28–34] . These porous urban flood models use a single isotropic porosity to account for the buildings in the cell, assuming
isotropic behavior. The reason for this assumption is that the shallow water equations with single porosity are derived from
the differential form of the classical shallow water equations using a representative elementary volume (REV) similar to
the derivation of the Darcy flow equation in groundwater flow modeling, e.g. [35] . The REV is by definition isotropic and
therefore only a single isotropic porosity can be derived for each cell (cf. [29] ), which leads to the loss of directionality and
hence may falsify the preferential flow paths. To the authors’ knowledge, two approaches have been developed to overcome
the loss of directionality and both have been developed for building treatment in urban flood models. The first approach has
been developed by Guinot [36] and introduces multiple porosities in each cell, which account for different directions and
storages. These porosities can be derived from the differential form of the shallow water equations without violating the
continuum model and REV assumption. The second approach has been introduced by Sanders et al. [37] . This approach ad-
ditionally assigns so-called areal or conveyance porosities to the cell edges, which introduce directionality to the equations.
If the differential form of shallow water equations is used, these areal porosities cannot be introduced without violating the
REV assumption. Therefore, the integral form of the shallow water equations is used, as it does not require the assumption
of an REV for the derivation. Yet, using the integral form of the shallow water equations means that only a finite volume
method can be utilized for the numerical solution [36,37] . Because these types of models are not isotropic anymore, they
are referred to as anisotropic porosity shallow water models.
While there is ongoing research at the University of Liege to incorporate depth-dependent porosities into an urban flood
model [38] , the porous shallow water models for building treatment generally do not allow full inundation of the buildings.
This is a valid assumption for urban flood modeling, however a porous shallow water model for generalized flow requires
partial as well as full inundation of unresolved topography. Therefore, this paper examines the possibility of extending
the equations derived in [37] to enable full inundation of the subgrid-scale unresolved topography to apply it to general
surface flow modeling. This leads to a novel formulation of the porosities and the interfacial pressure terms and a mutual
dependency between water elevation and porosity.
Finally, it should be noted that the shallow water equations with porosity cannot reproduce a high-resolution solution
exactly, because they cannot resolve local details of the flow. However, the anisotropic porosity model has been found to be
able to reproduce overall flow characteristics with satisfactory accuracy.
This paper is organized as follows: first the integral shallow water equations with anisotropic porosity are presented
in Section 2 ; then the numerical methods are discussed briefly; five computational examples are shown to demonstrate
model capabilities in Section 3 ; finally conclusions are given in Section 4 . In the following, the unresolved subgrid-scale
topography features such as microtopography in overland flow modeling or buildings in urban flood modeling are referred
to as unresolved solid structures or unresolved topography.
2. Governing equations
In this section, the integral shallow water equations with anisotropic porosity are derived for an arbitrary control volume.
As aforementioned, the numerical solution of the shallow water equations in integral form is only possible with the finite
volume method.
2.1. Anisotropic porosities
For the derivation of porosity, the phase function i inside a given control volume is introduced as:
i
(
x, y
)
=
½1 , if η(
x, y
)
> z
b
(
x, y
)
0 , else
(1)
where ηis the water elevation, z
b
is the bottom elevation and x, y are the horizontal Cartesian coordinates. Fig. 1 illus-
trates the water elevation ηand bottom elevation in a vertical section through a control volume. For illustration purposes,
i is evaluated in two points. If the bottom elevation exceeds the water elevation, i.e. dry case, the phase function is 0. If
the water elevation exceeds the bottom elevation, i.e. wet case, the phase function equals 1. Therefore, the phase function
indicates whether a certain point ( x, y ) in the control volume is wet or dry. Porosity is defined as the ratio of volume or
area of fluid to the whole volume or area of the control volume. Then, the volumetric porosity φis defined as:
φ= R
Äi
(
η−z
b
)
dÄ
R
Ä(
η−z
0
)
dÄ,(2)
2
Fig. 1. Definition of phase function i , water elevation η(blue), bottom elevation z
b (black) and zero datum z
0 (dashed) in a vertical section through a
control volume. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Fig. 2. Definition of control volume area Ä, control volume boundary ∂Äand path r in three dimensional view (left) top view and vertical section through
the cell edge marked with A–A (right). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this
article.)
and the areal porosity ψ of the boundary of the control volume is defined as:
ψ = H
∂Äi
(
η−z
b
)
dr
H
∂Ä(
η−z
0
)
dr
,(3)
where Ästands for the area of the control volume, ∂Ästands for the boundary of the control volume, z
0
is the zero datum
of the control volume (cf. Fig. 1 , dashed line) and r is the path along the boundary ∂Ä. The values of the porosities depend
on the zero datum z
0
. Here, the lowest bed elevation inside the control volume (denoted as minimum in Fig. 1 ) is chosen
as the zero datum. This means, that in the finite-volume method the zero datum will vary for each cell.
Fig. 2 shows an exemplary control volume with the definition of Äand ∂Ä. Fig. 2 (left) shows a three-dimensional
view of a partially inundated control volume, where blue color indicates the water column and gray color indicates bottom
topography. Fig. 2 (right) shows the top view of the control volume (top) and a vertical section through the boundary of
the control volume denoted with (A–A
′
) (bottom). Here, darker shades of gray indicate higher bottom elevation. Again, in
each point where the water inundates the bottom topography the phase function i = 1 and at the points where the bottom
topography elevation exceeds the water elevation i = 0 . Both elevations are calculated with the minimum bottom elevation
inside the control volume as the zero datum, which is marked in Fig. 2 (left). The volumetric porosity φis calculated with
3
Eq. 2 , and is the ratio of the volume of the fluid (blue columns) to the volume of the control volume. The volume of the
control volume is calculated by multiplying the elevation of the water column, i.e. the distance (A
′
–B
′
) in Fig. 2 (right,
bottom) with the total horizontal area of the cell Ä, shown in Fig. 2 (right, top). For example, in the case illustrated in
Fig. 2 , Ä= (A–A
′
)
2
. Similarly, the areal porosity ψ is calculated as the ratio of the vertical area of the fluid at the boundary
edge (colored blue in Fig. 2 (right, bottom)) to the vertical area of the boundary, described by the path (A–A
′
–B′
–B
′
) in Fig. 2
(right, bottom).
It can be shown that the constant porosities derived in [37] can be obtained by simplifying Eqs. 2 and 3 (cf. Appendix A ).
In contrast to these constant porosities, the porosities derived in this work are variable in time.
2.2. Integral-differential form of the shallow water equations with anisotropic porosity
The integral formulation of the shallow water equations can be obtained by applying the balance equation for mass and
momentum to a fixed Eulerian control volume under the assumption of hydrostatic pressure distribution [39] (pp. 47 ff.).
Prior to the integration, the conserved variables h, q
x and q
y are multiplied with the phase function i ( Eq. 1 ) to account for
the unresolved topography. Then, the temporal change of the vector of conserved variables q can be expressed as:
∂
∂t
Z
Ä
i q dÄ+
I
∂Ä
i Fn dr =
Z
Ä
s dÄ+
I
∂Ä∗
s
∗dr . (4)
Here, t is the time and s is the source vector. q and s are usually expressed as:
q =
"
h
q
x
q
y
#,s =
"i
r
s
b,x
+ s
f,x
s
b,y
+ s
f,y
#,(5)
where h = η−z
b
stands for water depth. q
x and q
y are the unit discharges in x - and y -direction, respectively. i
r is a mass
source term, e.g. rainfall intensity; s
b,x
, s
b,y
, s
f,x
, s
f,y
are the bed slope and the friction source terms in x - and y -direction,
respectively, and are calculated as:
s
b,x
= −gh
∂z
b
∂x
,s
b,y
= −gh
∂z
b
∂y
,(6)
s
f,x
= −c
f
q
x
pq
2
x
+ q
2
y
h
2
,s
f,y
= −c
f
q
y
pq
2
x
+ q
2
y
h
2
.(7)
The slope source terms account for variations in bottom and the friction source terms account for the bottom roughness. c
f
stands for the friction coefficient. F is the flux vector and can be expressed via f and g as:
Fn = f n
x
+ g n
y
, (8)
where n is the unit normal vector to the boundary; n
x and n
y are its components and f and g are the flux vectors in x - and
y -direction defined as:
f =
" q
x
uq
x
+ 0 . 5 gh
2
uq
y #
,g =
"q
y
v q
x
v q
y
+ 0 . 5 gh
2
#.(9)
Here, u and v are the depth-averaged velocity in x - and y -direction, respectively. g is the gravitational acceleration. s
∗is
the source vector accounting for fluid pressure along the interface between the fluid and solid ∂Ä∗. It results from the
macroscopic description, which does not differentiate between fluid and solid (cf. [40] , pp. 200–201).
In the limit of no structures to account for, the phase function i returns 1, the integral along ∂Ä∗vanishes and therefore
Eq. 4 converges to the classical two-dimensional shallow water equations, which can be found in, e.g. [39] (p. 47). In [11] ,
while properties of the differential form of the equations are discussed, it is argued that the equations may fail to give a
good approximation for very shallow flow, because some of the assumptions made for the derivation, e.g. a smooth free
surface, do not hold. Although the assumptions made in deriving the integral form are not violated during very shallow
flow, other statements made in [11] still apply and may lead to an inaccurate approximation. Namely, very shallow flow
with partially dry area is dominated by the effects of bottom irregularities which direct most of water laterally which
increases the flow path and the amount of dissipated energy [11] . If these bottom irregularities are only conceptually taken
into account by using the porosity and the interfacial pressure terms, the model will not be able to reproduce the correct
flow paths and may underestimate the dissipated energy. Unresolved topography which lies inside the computational cell
can only be accounted for with the volumetric porosity. This is a limitation of the model, because directionality is introduced
in the model in form of the areal porosities, to which the unresolved topography inside the cell cannot contribute. Hence,
structures which would have influenced the flow direction, e.g. roads and curbs, but lie completely inside the cell, will not
affect the flow direction. As a result, their impact on the flow may be underestimated by the model.
4
2.3. Storage and flux terms
The porosity terms in Eqs. 2 and 3 are used to express discrete forms of the integral terms containing the phase
function i .
The evaluation of the integral of i q in Eq. 4 is considered. In the following, volume-averaged variables will be used to
find a suitable approximation for this integral. The volume-averaged water elevation is calculated as:
¯η= R
ÄiηdÄ
R
ÄidÄ.(10)
The volume-averaged velocity is calculated as:
¯
v = R
Äih v dÄ
R
ÄihdÄ,(11)
The volume-averaged variables are constant within the control volume Ä. Applying Eq. 10 to Eq. 2 and using η−z
b
= h
leads to:
φ= R
Äi
(
η−z
b
)
dÄ
R
Ä(
η−z
0
)
dÄ= R
ÄihdÄ
R
Ä(
¯η−z
0
)
dÄ.(12)
As established above, ¯ηis constant inside the control volume. Hence, the expression (
¯η−z
0
) is also constant inside the
control volume and can be taken outside of the integration:
φ= R
ÄihdÄ
(
¯η−z
0
)
R
ÄdÄ= R
ÄihdÄ
(
¯η−z
0
)
Ä.(13)
Then, Eq. 13 can be rearranged to:
Z
Ä
ihdÄ=φ(
¯η−z
0
)
Ä, (14)
which corresponds to the evaluation of the integral of the mass storage (first entry of q ) in Eq. 4 . The momentum storage
in x -direction (second entry of q ) can be written by using q
x
= uh as:
Z
Ä
iq
x
dÄ=
Z
Ä
iuhdÄ. (15)
If the velocity u is approximated by the volume-averaged velocity, the equation becomes:
Z
Ä
iuhdÄ≈Z
Ä
i
¯
u hdÄ= ¯
u
Z
Ä
ihdÄ. (16)
Then, Eq. 13 can be used to write:
¯
u
Z
Ä
ihdÄ= φ¯
u
(
¯η−z
0
)
Ä. (17)
The same derivation can be applied in y -direction (third entry of q in Eq. 4 ) to get:
Z
Ä
iq
y
dÄ=
Z
Ä
i v hdÄ≈Z
Ä
i
¯
v hdÄ=
¯
v
Z
Ä
ihdÄ= φ¯
v
(
¯η−z
0
)
Ä. (18)
The integral of q in Eq. 4 can be replaced using Eqs. 13, 17 and 18 to write:
∂
∂t
φį
q +
I
∂Ä
i Fn d r =
Z
Ä
s d Ä+
I
∂Ä∗
s
∗dr, (19)
where ¯
q is the vector of volume-averaged variables:
¯
q =
"
(
¯η−z
0
)
¯
u
(
¯η−z
0
)
¯
v
(
¯η−z
0
)
#.(20)
The integral of i Fn in Eq. 4 can be evaluated by defining the area-averaged variables. Here, the area under consideration
( ∂Ä) is the boundary of the control volume. The closed curve integral of an arbitrary variable q can be splitted into integrals
along n segments:
I
∂Ä
qd r =
Z j+1
j
qd r +
Z j+2
j+1
qd r + · · · +
Z j
j+ n
qd r . (21)
This is illustrated in Fig. 3 , where the black line denotes ∂Äand the blue line denotes a piecewise linear approximation of it.
The approximation is intentionally crude for a better illustration. In theory, the integration can be carried out on the splitted
5
Fig. 3. Definition of path index k and vertex index j : top view of an arbitrary control volume; black color indicates the exact boundary, blue color indicates
the approximated boundary. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
parts of ∂Ä( Fig. 3 , black line), however in a finite volume method context the integration is carried out on piecewise linear
approximations of the boundary ( Fig. 3 , blue line). The area-averaged variables are calculated as:
ˆ
h = R
r
ihdr
R
ridr
,(22)
ˆ η= R
r
iηdr
R
ridr
,(23)
ˆ
v = R
r
ih v dr
R
rihdr
,(24)
r is the path between two points on ∂Ä, measured counter clockwise around ∂Ä, e.g. the path between point j and j + 1 in
Fig. 3 (marked with the index k + 1 ). The relationship between ∂Äand r is that ∂Äis the sum of all paths r .
ˆ
h is the area-
averaged water depth, ˆ ηis the area-averaged water elevation and ˆ
v =
¡ˆ
u ,
ˆ
v
¢is the area-averaged velocity vector. In a finite
volume method, variables would be averaged per cell edge, thus the area would be the edge under consideration. Therefore,
the term edge-averaged value is used interchangeably. To differentiate the area-averaged values from the volume-averaged
values, the area-averaged values are denoted with a circumflex (hat), e.g.
ˆ
h , and the volume averaged values are denoted
with a bar, e.g.
¯
h .
The flux term Fn in Eq. 4 is:
Fn =
q
x
n
x
+ q
y
n
y
¡uq
x
+ 0 . 5 gh
2
¢n
x
+ v q
x
n
y
uq
y
n
x
+
¡v q
y
+ 0 . 5 gh
2
¢n
y
.(25)
Eq. 3 can be rearranged to:
Z
r
i
(
η−z
b
)
d r = ψ
Z
r
(
η−z
0
)
d r . (26)
Further, applying the relation h = η−z
b
and Eq. 23 leads to:
Z
r
ihd r = ψ
Z
r¡ˆ η−z
0
¢d r = ψ
¡ˆ η−z
0
¢r . (27)
Eq. 24 in combination with Eq. 27 can be rearranged to:
Z
Ä
ih v d r =
ˆ
v
Z
Ä
ihd r = ψ
ˆ
v
¡ˆ η−z
0
¢r . (28)
This can be written in x - and y -direction as:
Z
r
ihudr = ψ
ˆ
u
¡ˆ η−z
0
¢r , (29)
6
Fig. 4. Definition of the interface ∂Ä∗(blue); gray blocks represent elements of microtopography. (For interpretation of the references to color in this
figure legend, the reader is referred to the web version of this article.)
and
Z
r
ih v dr = ψ
ˆ
v
¡ˆ η−z
0
¢r, (30)
respectively. Using q
x
= hu and q
y
= h v , the integral of the mass flux (first entry of Fn ) is approximated as:
Z
r
(
q
x
n
x
+ q
y
n
y
)
dr = ψ
ˆ
u
¡ˆ η−z
0
¢rn
x
+ ψ
ˆ
v
¡ˆ η−z
0
¢rn
y
. (31)
The momentum fluxes (second and third entries of Fn ) are approximated by using the area-averaged values
ˆ
h , ˆ
u and ˆ
v . In
x -direction this results in:
Z
r¡ih
ˆ
u
ˆ
u n
x
+ 0 . 5 igh
2
n
x
+ ih
ˆ
v
ˆ
u n
y
¢dr . (32)
The area-averaged values are taken outside of the integral:
ˆ
u
ˆ
u n
x
Z
r
ihd r +
Z
r
0 . 5 ig
ˆ
h hn
x
d r +
ˆ
v
ˆ
u n
y
Z
r
ihd r . (33)
Eq. 27 can be used to rewrite Eq. 32 :
ˆ
u
ˆ
u n
x
ψ
¡ˆ η−z
0
¢r + 0 . 5 g
ˆ
h n
x
ψ
¡ˆ η−z
0
¢r +
ˆ
v
ˆ
u n
y
ψ
¡ˆ η−z
0
¢r . (34)
The approximation of the momentum flux in y -direction is straight forward. Using Eq. 21 to replace the closed curve integral,
Eq. 19 is rewritten as:
∂
∂t
φį
q +
X
k
ψ
k
r
k
ˆ
F
k
n
k
=
Z
Ä
s dÄ+
I
∂Ä∗
s
∗dr, (35)
where k is the index of the path integral. The vector ˆ
F n is written as:
ˆ
F n =
ˆ
u
¡ˆ η−z
0
¢n
x
+
ˆ
v
¡ˆ η−z
0
¢n
y
ˆ
u
ˆ
u
¡ˆ η−z
0
¢n
x
+ 0 . 5 g
ˆ
h
¡ˆ η−z
0
¢n
x
+ ˆ
u
ˆ
v
¡ˆ η−z
0
¢n
y
ˆ
v
ˆ
u
¡ˆ η−z
0
¢n
x
+
ˆ
v
ˆ
v
¡ˆ η−z
0
¢n
y
+ 0 . 5 g
ˆ
h
¡ˆ η−z
0
¢n
y
.(36)
2.4. Solid–fluid interfacial pressure source term
∂Ä∗is the interface between fluid and solid, denoted with blue lines in Fig. 4 , where the top view of a square-shaped
control volume is given. The dashed black line shows the boundary of the control volume ( ∂Ä) and the gray blocks represent
single elements of simplified structures, e.g. buildings. Representing the unresolved fluid pressure at the interface ∂Ä∗, s
∗
consists of two components; the stationary component s
∗
st
which can be calculated if hydrostatic pressure distribution at
7
Fig. 5. Definition of p
∗and z
∗
b
; partially submerged control volume (left), fully submerged control volume (right): blue color indicates the water column.
(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
the interface is assumed and the non-stationary component s
∗
ns
which accounts for drag effects of the unresolved structures
[37] :
I
∂Ä∗
s
∗d r =
I
∂Ä∗
s
∗
st
d r +
Z
Ä
s
∗
ns
d Ä. (37)
While the stationary component s
∗
st
acts along the interface ∂Ä∗, the non-stationary component acts on the whole control
volume Ä.
In theory, the calculation of the stationary component s
∗
st
is straight-forward. Fig. 5 shows a vertical section through a
control volume and the two possible cases of submergence: partially submerged (left) and fully submerged (right). If these
cases are considered separately and hydrostatic pressure is assumed, the pressure of the fluid on the solid p
∗can be written
as:
p
∗(
x, y
)
=
(0 . 5 g
(
η−z
b
)
2 if η(
x, y
)
≤z
∗
b
,
0 . 5 g
¡¡η−z
∗
b
¢+
(
η−z
b
)
¢¡z
∗
b
−z
b
¢,else
(38)
where z
∗
b
is the bottom elevation of the microtopgraphy that the fluid pressure is acting on (cf. Fig. 5 ). If m =
¡m
x
, m
y
¢is
the unit normal vector along ∂Ä∗, which points inside the solid structure as illustrated in Fig. 4 , the stationary component
of the interfacial pressure source term can be written as:
s
∗
st
=
"0
p
∗m
x
p
∗m
y
#.(39)
The difficulty in the calculation of s
∗
st
is that the interface between solid and fluid ∂Ä∗is unknown because it is not resolved.
Therefore, the stationary term cannot be solved exactly and has to be approximated. One approach to estimate s
∗
st
can be
found in [37] .
The non-stationary component of Eq. 37 essentially accounts for drag which occurs during the flow through the un-
resolved structures, e.g. buildings or microtopography, as the fluid moves between the single elements of the structure.
Because it occurs at unresolved scales, the drag force cannot be calculated. In [37] , a generalized drag law is suggested to
approximate this term as:
s
∗
ns =
" 0
c
D
u | v |
c
D
v | v |#.(40)
Here, | v | is the Euclidian norm of the vector of velocities v =
(
u, v
)
and c
D
is a dimensionless drag coefficient. The determi-
nation of c
D
is challenging, often requires a calibration process and has not been fully understood yet. Several approaches
have been suggested to overcome this difficulty. In [37] , it is acknowledged that the drag effect may be estimated by an
increased roughness coefficient as demonstrated in [7] . In [11] , momentum correction terms are calculated which depend
8
on the volumetric porosity and a so-called effective water depth, which is the water volume per unit area. Also, differ-
ent methods with varying complexity for estimating c
D
have been presented in [28–30,37] . In [37] , a vegetative resistance
model as proposed in [41] is used to estimate c
D
. In this study, the drag force approach is used, because it is commonly
used and studied in literature, e.g. [28–30,37] . Here, the drag force approach of Nepf [41] is slightly modified to allow the
full submergence of the control volume. Then, c
D
is calculated as:
c
D
= 0 . 5 c
0
D
a ·min
¡h, z
∗
b
−z
b
¢.(41)
Here, a is the horizontally projected area of the elements of the solid structure per unit volume in one cell and c
0
D
is a bulk
drag coefficient accounting for the whole solid structure (cf. [41] ) and min is the minimum function. A similar approach is
given in [28] to account for inundated subgrid-scale structures. Both a and c
0
D
are not fully understood yet [37] , they depend
on the configuration of the solid structures as well as the shape of single elements, flow direction and several other factors
which have yet to be identified. Therefore, in this work the model is calibrated with the product c
0
D·a as a whole. Hence,
both a and c
0
D
lose their strict physical interpretations and become calibration parameters.
3. Numerical method and computational examples
The numerical solution of Eq. 35 can only be achieved with the finite volume method, as the equation does not contain
spatial differential expressions. Numerical studies of the authors have shown that a second order reconstruction of the
bottom elevation is necessary to obtain accurate results, especially in sloped domains (cf. [14] ). Further, a second-order
accurate scheme allows to compensate to some degree the loss of accuracy in the approximation due to coarse cells in
the anisotropic porosity model. Thus, for the following computational examples, the presented equations are solved with a
second-order monotonic upstream-centered scheme for conservation laws (MUSCL) presented in [42] being used for both
the anisotropic porosity model and the high-resolution model. It is acknowledged that the high-resolution model is suffering
more from the additional calculations per cell associated with the second order reconstruction process in comparison to the
anisotropic porosity model. A two-step explicit Runge–Kutta method is used to advance in time [43] . The numerical scheme
is implemented in the Hydroinformatics M odeling S ystem (hms), an in-house scientific programing framework [14] .
3.1. Calculation of porosities
Similar to [11] , it is suggested to calculate the porosities φand ψ with statistical properties of the unresolved subgrid-
scale features of the topography.
In a preprocessing step, the bottom elevation in each computational cell is sampled on a finer scale such that the discrete
cumulative distribution function (CDF) can be calculated individually in each cell. The CDF can then be used at the beginning
of each time step to evaluate how many of the samples are submerged by the water depth inside the cell. Basically, the CDF
is used as a structure to store the different bottom elevations mapped to the number of their occurrences.
For example, let the computational cell have a CDF based on 25 samples of bottom elevation inside the cell. It is assumed
that each sample corresponds to an equal area inside the cell. For sake of simplicity, let 10 of the samples have a bottom
elevation of 0, let 10 of the samples have a bottom elevation of 0.2 m and let 5 of the samples have a bottom elevation of
0.4 m. Then, the zero datum in the cell is defined as 0 and the water depth h in the cell corresponds to the water elevation
ηas h = η−0 = η. Further, assume the water depth is h = 0 . 1 m . The volume of the water inside the cell corresponds to
V
w
= 10 ·0 . 1 ·c, where c is the area of one sample. The total volume is V
t
= 25 ·0 . 1 ·c. Then, the volumetric porosity is
calculated as φ(h = 0 . 1) = V
w
/V
t
= 10 / 25 .
If the water depth rises to h = 0 . 3 m , the volume of water becomes V
w
= 10 ·0 . 3 ·c + 10 ·(0 . 3 −0 . 2) ·c, and the total
volume becomes V
t
= 25 ·0 . 3 ·c. Hence, the volumetric porosity is calculated as φ(h = 0 . 2) = V
w
/V
t
= 40 / 75 .
The same approach is applied to calculate areal porosities.
3.2. Error and speedup calculation
In the following, computational examples are presented to evaluate the capability of the equations. To the authors’ knowl-
edge, no analytical solutions for the shallow water equations with anisotropic porosity have been reported in literature. The
shallow water equations with isotropic porosity in [11] are compared with large-scale real case applications. The analytical
and semi-analytical solutions presented in [29] are valid for isotropic porosity only. In [37] , the anisotropic porosity model
results are compared with measurement data.
Therefore, in this work four examples are presented where the high-resolution shallow water model (HR) results are
considered to be the reference solution. In a final example, the anisotropic porosity model (AP) results are compared with
measurement data. The resolution of the HR model is always chosen such that further refinement does not change the
result. Turbulence and fluid viscosity are neglected in all test cases.
In order to assess the quality of the model results, the L
1
-norm, defined as:
L
1
=
1
N
N
X
i
| x
i
−˙
x
i
| ,(42)
9
Fig. 6. Dam break on dry bed with sine-wave shaped microtopography: initial conditions for η0
= 0 . 03 m (left); bottom elevation distribution inside one
AP model cell and the HR model discretization (dashed lines) (right).
is used, where N stands for the total number of solutions, x
i
is the reference solution, ˙
x
i
is the model solution and i is the
sample index.
The computational benefit of the anisotropic porosity model is quantified using the speedup, defined as:
ζ=
t
HR
t
AP
,(43)
whereby t
HR
and t
AP
are the wall-times of the HR model and the AP model, respectively.
3.3. One-dimensional dam break on dry bed with sine-wave shaped microtopography
In this computational example, a one-dimensional dam break on dry bed is simulated. The domain is 6 m long, 0.5 m
wide and the initial water elevation is defined as:
η(
x
)
=
½η0
, x ≤3 m
z
b
(
x
)
, x > 3 m
(44)
η0
stands for the initial water elevation and is varied from 0.025 m to 0.06 m for different simulation runs. The bottom
elevation of the domain is described with a sine-wave as:
z
b
(
x
)
= A sin
³2 π
λx +
π
2
´+ 0 . 01 . (45)
Here, λis the wavelength and A is the amplitude of the sine-wave. In this example they are set to λ= 0 . 05 m and A =
0 . 01 m . Fig. 6 (left) shows the initial conditions for η0
= 0 . 03 m . Only the section from 2.5 m < x < 3.5 m is plotted because
the small wavelength of the sine-wave makes it difficult to illustrate the bottom elevation over the whole domain. At the
outlet of the domain at x = 6 m , an open boundary forcing the water elevation gradient to zero is set. All other boundaries
are closed boundaries. Bottom roughness is accounted for with a Manning’s coefficient of n = 0 . 016 sm
−1 / 3
.
The reference solution is obtained by using a classical shallow water model with an element size of 1x = 0 . 01 m (HR
model). As shown in Fig. 6 (right), this resolution is sufficient to explicitly discretize the bottom elevation ( Eq. 45 ). In con-
trast, the model with anisotropic porosities (AP model) uses a mesh with element size of 1x = 0 . 1 m . Fig. 6 (right) shows
exactly one computational cell of the AP model and the bottom topography inside it. The resolution of the AP model’s mesh
is not sufficient to explicitly discretize the sine-wave, therefore the bottom elevation is described by z
b
(
x
)
= 0 m . A classical
shallow water model (SWE model) with the same resolution as the AP model is used to illustrate the effect of the AP model.
Good agreement between the HR model and the AP model is achieved for c
0
D
·a = 10 m
2 as shown in Figs. 7 and 8 . In [37] ,
c
0
D
= 2 is recommended but values up to c
0
D
= 6 have been reported (M. Bruwer, personal communication, 24 March 2015)
which shows that this value has an uncertainty. The results for water elevation, velocity and unit discharge are plotted on
the left side in Fig. 7 (top, middle and bottom, respectively). The fluctuations in the HR model solution are due to the sine-
wave shaped microtopography as the water accelerates when flowing down the sine-wave and decelerates when climbing
up the crests of the sine-wave. The SWE model shows poor agreement in all cases. The AP model captures the advance
of the front correctly and the obstructive effects of the microtopography could be reproduced well ( Fig. 7 (top left)). The
velocity is underestimated between 0 < x < 1 m and slightly overestimated between 1 m < x < 4 m ( Fig. 7 (middle left)).
The unit discharge behaves similar as the velocity ( Fig. 7 (bottom left)). Water elevation, velocity and unit discharge are all
captured well. On the right side in Fig. 7 , the sensitivity of c
0
D·a is illustrated. The product c
0
D·a is varied from 0 to 500 m
2
.
As c
0
D
·a increases, the roughness of the model increases. The AP model is sensitive with regard to c
0
D·a until a critical
value of about c
0
D·a = 500 m
2 is reached. It was observed that for c
0
D·a > 500 m
2 this parameter is not very sensitive. This
10
Fig. 7. Dam break on dry bed with sine-wave shaped microtopography: comparison of model results at t = 4 s with η0
= 0 . 06 m (the anisotropic porosity
model (AP), the high-resolution reference solution (HR) and a coarse grid classical shallow water model (SWE)) (left), sensitivity study of c
0
D·a (denoted
as c) (right).
is because after reaching a certain value, friction is artificially limited in the numerical scheme to avoid velocities to change
direction. For details on this friction treatment, the reader is referred to [44] .
The initial water elevation η0
is varied to 0.025 m and 0.03 m to study the influence of the water elevation. Fig. 8 (left)
and (right) show that the solution is enhanced by the AP model for different initial water elevations. In both cases, the
overestimation of the velocity and the discharge is higher than for η0
= 0 . 06 m . It is noted, that the drag coefficient c
0
D
·a =
10 m
2 is kept constant for these simulations. Case specific calibration might further enhance the solution. The L
1
-error for
the presented cases is summarized in Tables 1 –3 for water elevation, velocity and discharge. For all variables, the L
1
-error
of the SWE model is about one order of magnitude higher than the AP model error. The computation with the AP model
was carried out approximately 10 0 0 times faster than with the HR model.
3.4. Two-dimensional dam break across a porosity discontinuity
The following example simulates a two-dimensional dam break across a porosity discontinuity on a 100 m ×10 m do-
main. This example was initially introduced in [29] as a one-dimensional benchmark for the shallow water model with
11
Fig. 8. Dam break on dry bed with sine-wave shaped microtopography: comparison of model results at t = 4 s with η0
= 0 . 03 m (left) and
η0
= 0 . 025 m
(right) (the anisotropic porosity model (AP), the high-resolution reference solution (HR) and a coarse grid classical shallow water model (SWE)).
Table 1
Dam break on dry bed with sine-
wave shaped microtopography: L
1
-error
for water elevation as calculated by the
anisotropic porosity model (AP) and the
coarse shallow water model (SWE).
η(m) L
1, AP
(
η) (m) L
1, SWE
(
η) (m)
0.025 3 . 8 ·10
−4 28 ·10
−4
0.03 6 . 4 ·10
−4 34 ·10
−4
0.06 15 ·10
−4 76 ·10
−4
12
Table 2
Dam break on dry bed with sine-wave
shaped microtopography: L
1
-error for veloc-
ity as calculated by the anisotropic poros-
ity model (AP) and the coarse shallow water
model (SWE).
η(m) L
1, AP
( v ) (m/s) L
1, SWE
( v ) (m/s)
0.025 0 . 8 ·10
−2 5 . 7 ·10
−2
0.03 1 . 6 ·10
−2 7 ·10
−2
0.06 1 . 2 ·10
−2 13 . 8 ·10
−2
Table 3
Dam break on dry bed with sine-wave shaped
microtopography: L
1
-error for unit discharge
as calculated by the anisotropic porosity
model (AP) and the coarse shallow water
model (SWE).
η(m) L
1, AP
( q ) (m
2
/s) L
1, SWE
( q ) (m
2
/s)
0.025 1 . 0 ·10
−4 8 . 6 ·10
−4
0.03 2 . 6 ·10
−4 12 ·10
−4
0.06 4 . 9 ·10
−4 44 ·10
−4
Fig. 9. Dam break across a porosity discontinuity: top view on the computational domain (left); top view on an exemplary computational cell of the AP
model (right), black color indicates the location of obstacles.
isotropic porosity and a quasi-analytical solution for the one-dimensional case was derived. This solution is not valid for
the two-dimensional case, therefore the results of the anisotropic porosity model (AP) is compared with a high-resolution
shallow water model (HR). The computational domain is illustrated in Fig. 9 (left). The discontinuity of the porosity as well
as the discontinuity of the water elevation is located at x = 50 m :
η0
(
x, y
)
=
½2 m , x ≤50 m
1 m , x > 50 m
φ0
(
x, y
)
=
½1 , x ≤50 m
0 . 8 , x > 50 m
(46)
At the outlet x = 100 m , an open boundary forcing the water elevation gradient to zero is set. All other boundaries are
closed wall boundary conditions. The porosity jump is constructed via randomly generated obstacles which are explicitly
discretized in the HR model and are taken into account by the porosities in the AP model. All obstacles are square shaped
with an edge length of 0.1 m and with infinitive vertical height and are spatially distributed according to a random uniform
distribution such that each cell of the AP model has a volumetric porosity of φ= 0 . 8 for x > 50 m as illustrated in Fig. 9
(right) for one exemplary cell. During the whole simulation, the obstacles are never fully inundated, which means that
the volumetric porosity stays constant in each cell. The simulation is run for t = 4 s . The HR model is calculated on a grid
with square-shaped elements with an edge length of 0.02 m. The AP model uses a computational grid with square-shaped
elements with an edge length of 0.5 m. c
0
D
·a = 10 m
2
is chosen for the AP model. Bottom roughness in both models is taken
into account by a Manning’s coefficient of n = 0 . 016 sm
−1 / 3
.
Results for water elevation and unit discharge at different longitudinal sections at t = 4 s are plotted in Fig. 10 (left) and
(right), respectively. L
1
-errors for water elevation and unit discharge at the sections are given in Table 4 . The AP model
results show good agreement with the reference solution calculated by the HR model. After the dam break at x = 50 m , the
rarefaction wave traveling in upstream direction as well as the shock wave traveling in flow direction are captured well,
13
Fig. 10. Dam break across a porosity discontinuity: water elevation (left) and unit discharge (right) at t = 4 s for different longitudinal sections for the
anisotropic porosity model (AP) and the high-resolution reference solution (HR).
although at about x = 30 m the water elevation is underestimated in all sections. The fluctuation of the water elevation,
which results from the superposition of waves due to the obstacles, calculated by the HR model cannot be reproduced by
the AP model. The unit discharge is captured very well by the AP model ( Fig. 10 (right)). The discretized obstacles in the
HR model narrow the cross section available to flow and lead to a high localized flow velocity and therefore a high unit
discharge. This cannot be reproduced by the AP model. As pointed out in [29] , this is not a failure of the AP model, but is a
consequence of the macroscopic modeling using the porosity concept. The AP model results were computed roughly 30 0 0
times faster than the HR model results.
3.5. Two-dimensional dam break on dry bed with random microtopography
This example considers a two-dimensional dam break on dry bed with random microtopography. The domain spans 6 m
in x -direction and 3 m in y -direction. The water elevation is defined as:
η(
x
)
=
½0 . 03 m , x ≤3 m
z
b
(
x, y
)
, x > 3 m
(47)
14
Table 4
Dam break on dry bed across a poros-
ity discontinuity: L
1
-error for water
elevation and unit discharge.
y (m) L
1
(
η) (m) L
1
( q ) (m
2
/s)
0.525 8 . 1 ·10
−3 5 . 3 ·10
−2
1.455 7 ·10
−3 3 . 9 ·10
−2
2.25 7 . 4 ·10
−3 5 . 3 ·10
−2
Fig. 11. Dam break on bed with random microtopography: initial conditions (left); microtopography in the domain (right).
Table 5
Dam break with random microto-
pography:
L
1
-error for water eleva-
tion and velocity for y = 0 . 525 m .
t (s) L
1
(
η) (m) L
1
( v ) (m/s)
0.4 1 . 8 ·10
−4 2 ·10
−3
1 2 . 5 ·10
−4 3 . 6 ·10
−3
1.8 6 ·10
−4 9 . 1 ·10
−3
Table 6
Dam break with random microto-
pography:
L
1
-error for water eleva-
tion and velocity for y = 0 . 525 m .
t (s) L
1
(
η) (m) L
1
( v ) (m/s)
0.4 2 . 7 ·10
−4 4 ·10
−3
1 2 . 5 ·10
−4 5 ·10
−3
1.8 2 ·10
−4 6 ·10
−3
The microtopography is generated as square-shaped deviations with an edge length of 0.05 m, and their amplitudes z
b,mic
are
distributed between 0 and 0.02 m according to a Gaussian distribution function as illustrated in Fig. 11 (right). All boundaries
are closed except at the right side of the domain ( x = 6 m ), where an open boundary condition as in previous the example
is applied. A reference solution is computed with a shallow water model on a 0.01 m ×0.01 m grid (HR). The anisotropic
porosity model uses square grid cells with an edge length of 0.1 m (AP). The bottom friction is expressed via a Manning
coefficient n = 0 . 016 sm
−1 / 3
. The drag force of the AP model is estimated with the product c
0
D·a = 10 m
2
. The simulation
runs for t = 2 s .
Results for water depth at different sections through different y -values are plotted in Fig. 12 . Here, dry cells are not
plotted for the HR model. The L
1
-errors for water elevation and velocity at different times are given in Tables 5 and 6 . The
AP model shows very good agreement with the HR model. The shock is captured with satisfactory accuracy at all times,
however local details of the water elevation variation such as small scale fluctuations due to the microtopography cannot
be captured.
Velocity profiles through the same sections as in Fig. 12 are plotted in Fig. 13 . Although the maximum values of the
velocity profiles are not reproduced by the AP model, overall good agreement between the HR model and the AP model is
observed.
15
Fig. 12. Dam break on bed with random microtopography: water elevations at y = 0 . 525 m (left) and y = 2 . 245 m (right) at different times for the
anisotropic porosity model (AP) and the high-resolution reference solution (HR); the high-resolution bottom topography is plotted at the very top of
each column for illustration purposes.
16
Fig. 13. Dam break on bed with random microtopography: flow velocities at y = 0 . 525 m (left) and y = 2 . 245 m (right) at different times for the anisotropic
porosity model (AP) and the high-resolution reference solution (HR).
Fig. 14 shows a top view on the water elevation distribution in the domain at the same time steps as in Fig. 12 for the HR
model (left) and the AP model (right). The HR model resolves the microtopography explicitly and as the water elevation is
calculated as η= h + z
b
, in the dry part of the domain, the water elevation equals the bottom elevation. It is observed that
the overall characteristics of the advancing front and the rarefaction wave moving upstream are captured well by the AP
model. However, the spatial distribution of the AP model results have low accuracy, as they suffer from numerical diffusion
due to coarse grids as well as the lack of information on small scale bottom elevation variations. The results of the AP model
are calculated approximately 10 0 0 times faster than the HR model.
3.6. Rainfall-runoff on an inclined plane with random microtopography
Rainfall-runoff is heavily influenced by the microtopography of the domain [45] . In this example, the surface runoff on
a 6 m ×3 m inclined plane with a slope of 0.02 and a Manning’s coefficient of n = 0 . 016 sm
−1 / 3 is simulated. The bottom
elevation for the high-resolution model (HR) is calculated as:
z
b
(
x, y
)
= 1 −0 . 02 ·x + z
b,mic
(
x, y
)
. (48)
17
Fig. 14. Dam break on bed with random microtopography: water elevations at different time steps for the high-resolution reference solution (HR) (left)
and the anisotropic porosity model (AP) (right).
Fig. 15. Rainfall-runoff on an inclined plane with random microtopography: side view of the computational domain without microtopography (left); top
view of the position and spatial distribution of microtopography (right top); top view of the positions of evaluation points (right bottom).
Here, z
b, mic
is the amplitude of the microtopography, which is generated as square blocks with an edge length of 0.02 m
and a vertical amplitude varying between 0 and 0.003 m according to a Gaussian distribution function. The microtopogra-
phy is applied only inside a rectangular area spanning from (2.25 m, 0.75 m) to (3.75 m, 2.25 m) whereby the first pair of
coordinates denotes the bottom left corner and the latter pair denotes the top right corner of the rectangle as illustrated in
Fig. 15 (right top). For the anisotropic porosity model (AP), the microtopography is not explicitly discretized and the bottom
elevation is calculated as:
z
b
(
x, y
)
= 1 −0 . 02 ·x . (49)
The domain without microtopography is illustrated in Fig. 15 (left). Rainfall is imposed for 100 s with the intensity being
varied from i
r
= 1 ·10
−5
m / s to i
r
= 1 ·10
−3
m / s for different simulation runs. The boundary at the outlet is an open bound-
ary, all other boundaries are closed. The HR model uses a square grid with an element size 1x = 0 . 02 m , the AP model uses
a square grid with an element size of 1x = 0 . 1 m . A calibration resulted in c
0
D·a = 0 , i.e. no drag force influence.
18
Fig. 16. Rainfall-runoff on an inclined plane with random microtopography: comparison of normalized discharges (left) and water depths (right) at the
outlet computed by the HR model and the AP model for different rainfall intensities.
The normalized discharges at the outlet of the domain ( x = 6 m ) are compared for the different rainfall intensities in
Fig. 16 . The normalized discharge is calculated as the ratio of the model discharge ( Q
model
) to the rainfall discharge, whereby
the rainfall discharge for a l ×w rectangular domain is calculated as [46] :
Q
rain
= l ·w ·i
r (50)
The comparison shows that the influence of the microtopography is overestimated by the AP model. In the early time of the
simulation, both hydrographs agree well but when the wave which is influenced by the microtopography reaches the outlet
the hydrographs start to deviate.
For i
r
= 1 ·10
−5
m / s the AP model does not reach its concentration time in 100 s. The agreement at the late stages of
the simulation (after t = 80 s ) is less good. This suggests that in the AP model, the influence of the microtopography is
overestimated in these test cases and thus the water is artificially held back and does not reach the outlet. This argument
is supported by the fact that the agreement gets better for increasing rainfall intensity, cf. e.g. the hydrograph of i
r
= 1 ·
10
−3
m / s . As the intensity increases, the influence of the microtopography on the flow decreases. For i
r
= 1 ·10
−4
m / s the
hydrograph of the AP model rises a little bit too slow and for i
r
= 1 ·10
−3
m / s both hydrographs agree well. The water
depths behave similarly. The results of the AP model are on average computed 550 times faster than the reference solution.
19
Table 7
Rainfall-runoff on an inclined plane with random
microtopography: scaled L
1
-error for unit discharge
and water elevation at the outlet, errors are scaled
by division by the corresponding rainfall intensity.
i (m/s) L
1
( q ) ((m
2
/s)/(m/s)) L
1
( h ) (m/(m/s))
10−31 . 4 ·10
−2 3 . 8 ·10
−2
10−410 ·10
−2 40 ·10
−2
10−528 ·10
−2 280 ·10
−2
Table 8
Rainfall-runoff on an inclined
plane with random microtopog-
raphy: scaled L
1
-error for unit
discharge at the gauges for i =
10
−4
m / s , errors are scaled by
division by the rainfall intensity.
Point L
1
( q ) ((m
2
/s)/(m/s))
1 2 ·10
−1
2 1 . 6 ·10
−1
3 10 ·10
−1
4 3 . 1 ·10
−1
5 0 . 13 ·10
−1
Table 9
Rainfall-runoff on an inclined
plane with random microtopog-
raphy: scaled L
1
-error for unit
discharge at the gauges for i =
10
−5
m / s , errors are scaled by
division by the rainfall intensity.
Point L
1
( q ) ((m
2
/s)/(m/s))
1 5 . 1 ·10
−1
2 7 . 2 ·10
−1
3 30 ·10
−1
4 4 ·10
−1
5 1 . 2 ·10
−1
The L
1
-errors for different intensities are given in Table 7 . Here, the L
1
-error is divided by the corresponding intensity
for better comparison of the cases.
For i
r
= 1 ·10
−4
m / s , model results at different points are compared (cf. Fig. 17 , right bottom). Fig. 17 also shows a com-
parison of normalized discharges at these evaluation points. Good agreement between the discharges is observed at the
points 1, 2 and 5. However, especially at points 1 and 2 is a temporal delay in the hydrograph of the AP model which
again comes from the overestimation of the influence of microtopography. Point 3, which is located inside the area with
microtopography, shows the worst agreement which might result from the aforementioned overestimation as well as the
macroscopic approach of the AP model which is not expected to reproduce local flow processes. At point 4 the discharge is
overshot by the AP model.
Model results for the same points are compared in Fig. 18 for i
r
= 1 ·10
−5
m / s . Here, it is seen that the agreement at
the points where the flow is influenced by the microtopography, namely points 2, 3 and 4, gets worse for lower intensities.
Especially at point 3, which is located inside the area with microtopography, the AP model returns a discharge which is 3
times higher than the HR model discharge and is temporally delayed.
The L
1
-errors at the different points are given in Tables 8 and 9 for i = 10
−4
m / s and i = 10
−3
m / s , respectively. Again,
the L
1
-errors are divided by the corresponding intensity.
Fig. 19 shows temporal snapshots of the discharge distribution in the domain at t = 15 s , t = 20 s and t = 50 s for both the
HR model (left) and the AP model (right). The resolution of the AP model is much coarser than the HR model and therefore
local details cannot be resolved as good as in the HR model but general properties of the flow field are reproduced. At
t = 20 s the overestimation of microtopography can be seen very clearly, as the flow calculated by the HR model ( Fig. 19 ,
left middle) has already reached the right border of the microtopography area while the flow calculated by the AP model
( Fig. 19 , right middle) has only reached the middle of the microtopography area. The discharge of the AP model is higher
20
Fig. 17. Rainfall-runoff on an inclined plane with random microtopography: comparison of normalized discharges computed by the HR model and the AP
model for i
r
= 1 ·10
−4
m / s at different evaluation points (plotted in the right bottom).
than of the HR model, however the porosity slows down the front of the AP model flow. This can also be observed in Fig. 17
(left middle), where the discharge at point 3 is delayed and overestimated by the AP model. At t = 50 s the flow fields
reasonably resemble ( Fig. 19 , bottom).
3.7. Dam-break flow through an idealized city
In this computational example, results of a dam-break experiment conducted at the Université catholique de Louvain,
Belgium, [47] are numerically reproduced.
3.7.1. Domain description, initial and boundary conditions
The domain is a 35.8 m long and 3.6 m wide channel with horizontal bed. The idealized city consists of 5 ×5 buildings,
each of them being a square block with an edge length of 0.30 m. The distance between the blocks is 0.10 m. The dam-
break is constructed by opening a 1 m wide gate, which initially separates part of the channel with water ponding at 0.40 m
from the rest of the channel (reservoir), where a very thin layer of 0.011 m water due to imperfect tightness of the gate
21
Fig. 18. Rainfall-runoff on an inclined plane with random microtopography: comparison of normalized discharges computed by the HR model and the AP
model for i
r
= 1 ·10
−5
m / s at different evaluation points (plotted in the right bottom of Fig. 17 ) .
is reported. For further details on the experimental setup and employed measurement techniques, the reader is referred to
[47] . The domain is illustrated in Fig. 20 , where the reservoir is colored in gray.
For the numerical model, only the reservoir and the first 16 m of the channel is discretized for computational efficiency. In
preliminary studies it had been observed that for the duration of the simulations, t = 15 . 5 s , the shock wave does not travel
further than this length. The downstream boundary is an open boundary and all other boundaries are closed boundaries.
The HR model uses a triangular mesh with three different cell sizes: the inside of the reservoir is discretized with cells
with a characteristic length of l
c, 1
= 0 . 3 m . The area inside the channel which is sufficiently far away from the building blocks
is discretized with a characteristic length of l
c, 2
= 0 . 1 m . The space between the buildings is discretized with a characteristic
length of l
c, 3
= 0 . 01 m . The buildings are represented as holes in the mesh, which is a method commonly used in urban
flood modeling [48] . With the value chosen for l
c ,3
, the space between two buildings is discretized with about 10 cells and
the total cell number is 95975. The AP model uses square-shaped cells with an edge length of 0.25 m in the whole domain,
which results in 1272 cells in total.
Experimental data is available at 64 measurement gauges distributed inside the channel [47] . The positions of these
gauges are given in Fig. 21 . Errors are calculated for all gauges. In the discussion, results are plotted only for 4 gauges,
namely gauges 1, 18, 44 and 55, to avoid too many figures. These gauges are pointed out in Fig. 21 .
The roughness of the channel has been estimated in [47] with a Manning’s coefficient of n = 0 . 01 sm
1 / 3
. This value is
used for both the HR and the AP model.
22
Fig. 19. Rainfall-runoff on an inclined plane with random microtopography: comparison of snap shots of unit discharges computed by the HR model (left)
and the AP model (right) at different time steps.
3.7.2. Model calibration and run time
The AP model is calibrated with the value a ·c
0
D
. The best results in this specific case were obtained by completely
neglecting the drag force, i.e. a ·c
0
D
= 0 . No calibration is carried out for the HR model. The HR model simulation takes about
30 0 0 s wall-clock time to finish. The AP model requires about 4 s wall-clock time. Consequently, the speedup is calculated
as 750 (cf. 10 ).
3.7.3. Discussion
The HR model makes overall a good prediction of the water depth at the evaluated gauges. In Fig. 22 , the water depth
calculated by the HR model at the aforementioned gauges is plotted together with the measured water depth. Overall, the
numerical results approximate the experimental results very well. The arrival time of the wave is predicted correctly at
all gauges. Larger deviations between the results occur at the later stages of the simulation. At gauge 18, which is located
between the buildings, the wave reflections from the walls of the buildings superpose and create several peaks between
t = 3 . 5 s and t = 6 . 5 s in the HR model results which were not observed during the experiment. Further, at gauge 1, which is
at the upper right corner of the building block, the water depth is underestimated by the HR model. This might be because
of the hydraulic jump observed at the impact section which leads to increased local water levels which are not reproduced
23
Fig. 20. Dam break flow through an idealized city: computational domain and initial conditions.
Fig. 21. Dam break flow through an idealized city: location of the gauges, area of building array is marked with dashed line.
Table 10
Summary of speedups obtained in all simulations, n : number of cells,
1x : edge
length, HR: high-resolution model, AP: anisotropic porosity model.
Case 1x
HR
(m) n
HR 1x
AP
(m) n
HR n
HR
/ n
AP Speedup
3.3 0 . 01 30 , 0 0 0 0 . 1 300 100 10 0 0
3.4 0 . 02 2 , 500 , 000 0 . 5 40 0 0 625 30 0 0
3.5 0 . 01 180 , 0 0 0 0 . 1 1800 100 10 0 0
3.6 0 . 02 45 , 0 0 0 0 . 1 1800 25 550
3.7 0 . 01 −0 . 3 95 , 975 0 . 25 1272 75.4 750
by the HR model. These points might raise the question, whether a turbulence model should be used, however Soares-Frazão
and Zech [47] report that adding turbulence to the numerical model leads to a worse agreement between numerical and
experimental results.
The anisotropic porosity model (AP) shows good agreement with the HR model results, although the results of the AP
model are smoother and more diffused than the HR model results. In Fig. 22 , AP model results for water depth are plotted
for the four gauges as well. Gauge 1 and gauge 18 show very good agreement, while the arrival time of the wave at gauge
44 is delayed. Gauge 55, located in the front of the building block, shows the worst agreement of the four. Here, the AP
model overshoots the HR model. The peak at around t = 4s is not reproduced. Overall, the general properties of the AP
model results, i.e. the lack of local and spatial fluctuations, agree with the observations in [49] . In general, the AP model
error manifests itself in excessive damping of the results.
L
1
-errors of both models in regard to the experimental data are calculated as the average L
1
-error of all 64 gauges. The
HR model has a L
1
-error of 0.02 m, the AP model has a L
1
-error of 0.07 m.
24
Fig. 22. Dam break flow through an idealized city: discharges calculated by the anisotropic porosity model (AP), high-resolution model (HR) and the
measurement data at gauges 1, 18, 44 and 55.
4. Conclusions
An integral formulation of the two-dimensional shallow water equations with anisotropic porosity for flow over partially
and fully inundated topography was derived. A novel formulation for the porosities was proposed and an approximation for
the storage and flux terms was presented. The porosities are dependent on the water elevation in the cell. This relationship
can be approximated by calculating a cumulative distribution function for the unresolved bottom elevation and evaluating
it at the water elevation. Due to the macroscopic point of view, additional terms appear in the governing equations. Suit-
able approximations for these terms have been referred to. The non-stationary term was approximated with a drag force
approach. The integral formulation of the equations can only be solved by the finite volume method. A second order MUSCL
scheme was used to solve the equations with a two-step explicit Runge–Kutta method for time stepping.
Five computational examples, ranging from simple academic benchmarks to nearly ‘real case’ laboratory experiments
were shown to demonstrate the capabilities and limitations of the new approach. Due to the lack of analytical solutions
a high-resolution shallow water model was used to calculate reference solutions. In the last test case, experimental data
was used for model evaluation. The shallow water model with porosity showed overall good agreement with the reference
solutions. The aforementioned drag term was used to calibrate the model and a sensitivity study regarding this term was
carried out. Except in the last test case, good results are obtained with c
0
D
·a = 10 . However, further studies to investigate
the drag force coefficient values and the possibility to represent the drag effect with increased friction are required.
As bottom slope increases, the accuracy of the anisotropic shallow water model decreases. Experimental studies show
that a large bed slope reduces the effect of microtopography [50] and the presented model seems to underestimate the
reduction.
A challenge in practical applications is the isolation of the part of topography modeled as porosity from the global to-
pography. Usually, the global topography is defined as the roughness of the surface of the earth and represented by the cell
value. The unresolved topography is thought about as subgrid-scale deviations from this value which creates heterogeneity
inside the cell. The issue of identifying these deviations has been researched in the context of isolating microtopography
and different methods have been proposed in the literature, e.g. [51] . However, finding suitable methods to correctly isolate
the part of topography to be modeled as porosity remains an open issue, which seems to be the main limitation of applying
the proposed model to ‘real world cases’.
Local details of the flow could not be exactly reproduced by the anisotropic porosity model, because the concept of
porosity as a statistical property of the topography is not expected to reproduce processes at this scale [29] .
The novel anisotropic porosity was found to be a good balance between computational time and accuracy. Table 10 gives
an overview of the speedups in the simulations in dependency of cell size 1x and cell number n . The ratio of cell numbers
( n
HR
/ n
AP
) is identified as the main factor of the speedup. In the presented computational examples, the anisotropic porosity
25
model provided a computational benefit around three orders of magnitude, depending on the ratio of the cell numbers, i.e.
the difference in cell size.
Acknowledgment
The research fellowship granted by the Alexander von Humboldt Foundation is gratefully acknowledged by Dongfang
Liang. The authors thank Prof. Andrea Defina, Università degli Studi di Padova, Italy, for the insightful personal communica-
tion, and the Local Organizing Committee of the Advances in Numerical Modelling of Hydrodynamics Workshop, Sheffield,
UK, guest editor Dr. Daniel Caviedes-Voullième and the anonymous reviewers for their valuable comments.
Appendix A. Derivation of porosities by Sanders et al. [37]
It can be shown that the definitions of porosity in [37] can be considered as a special case of Eqs. 2 and 3 , where
submergence of microtopography is not allowed. If microtopography is not allowed to be submerged, the wet fraction of
the control volume will remain constant. Further, if the wet fraction of the control volume is considered to have the same
constant bed elevation z
b
= z
0
, Eq. 2 can be simplified to:
φ= R
Äi
(
η−z
0
)
dÄ
R
Ä(
η−z
0
)
dÄ=
(
η−z
0
)
R
ÄidÄ
(
η−z
0
)
R
ÄdÄ=
1
ÄZÄ
idÄ. (A.1)
If the same assumptions are made for the boundary of the control volume, Eq. 3 simplifies to:
ψ = H
∂Äi
(
η−z
0
)
dr
H
∂Ä(
η−z
0
)
dr
=
1
∂ÄI∂Ä
idr. (A.2)
Eqs. A.1 and A.2 are the porosities introduced in [37] for building treatment. The assumptions made are valid for describing
the effects of buildings, which are unlikely to become submerged by the flood wave. If the porosities are used to describe
the effects of microtopography, Eqs. 2 and 3 have to be used. A significant difference between Eqs. 2 and 3 (with inundation)
and Eqs. A.1 and A.2 (without inundation) is that the porosities that allow inundation are dependent on the water elevation
η, which is variable in time. If the water elevation increases, the porosities increase. Therefore, the porosities are functions
of time if the terrain variation within a control volume is considered and inundation is allowed, while without inundation
the wet fraction of the control volume remains constant and thus, the porosities are constant in time.
References
[1] I. Özgen , S. Seemann , A.L. Candeias , H. Koch , F. Simons , R. Hinkelmann , Simulation of hydraulic interaction between Icó-Mandantes bay and São Fran-
cisco river, Brazil, in: G. Gunkel, J.A .A . Silva, M.d. C. Sobral (Eds.), Sustainable Management of Water and Land in Semiarid Areas, Editora Universitária
UFPE, Recife, Brazil, 2013, pp. 28–38 .
[2] G. Kesserwani, Q. Liang, RKDG2 shallow-water solver on non-uniform grids with local time steps: Application to 1D and 2D hydrodynamics, Appl.
Math. Model. 39 (3-4) (2015) 1317–1340, doi:
10.1016/j.apm.2014.08.009 .
[3] C. Lange , M. Schneider , M. Mutz , M. Haustein , M. Halle , M. Seidel , H. Sieker , C. Wolter , R. Hinkelmann , Model-based design for restoration of a small
urban river, J. Hydroenviron. Res. (JHER) 9 (2) (2015) 226–236 .
[4] G. Petaccia, L. Natale, F. Savi, M. Velickovic, Y. Zech, S. Soares-Frazão, Flood wave propagation in steep mountain rivers, J. Hydroinf. 15 (1) (2013) 120,
doi: 10.2166/hydro.2012.122 .
[5] C. Zoppou, S. Roberts, Numerical solution of the two-dimensional unsteady dam break, Appl. Math. Model. 24 (7) (20 0 0) 457–475, doi: 10.1016/
S0307-904X(99)0 0 056-6 .
[6] Y. Liu, J. Zhou, L. Song, Q. Zou, L. Liao, Y. Wang, Numerical modelling of free-surface shallow flows over irregular topography with complex geometry,
Appl. Math. Model. 37 (23) (2013) 9482–9498, doi: 10.1016/j.apm.2013.05.001 .
[7] D. Liang, R.A. Falconer, B. Lin, Coupling surface and subsurface flows in a depth averaged flood wave model, J. Hydrol. 337 (1-2) (2007) 147–158,
doi: 10.1016/j.jhydrol.2007.01.045 .
[8] D. Liang, X. Wang, R.A. Falconer, B.N. Bockelmann-Evans, Solving the depth-integrated solute transport equation with a TVD-MacCormack scheme,
Environ. Model. Softw. 25 (12) (2010) 1619–1629, doi: 10.1016/j.envsoft.2010.06.008 .
[9] Q. Liang, Flood simulation using a well-balanced shallow flow model, J. Hydraul. Eng. 136 (9) (2010) 669–675, doi: 10.1061/(ASCE)HY.1943-7900.
0 0 0 0219 .
[10] G. Petaccia, S. Soares-Frazão, F. Savi, L. Natale, Y. Zech, Simplified versus detailed two-dimensional approaches to transient flow modeling in urban
areas, J. Hydraul. Eng. 136 (2010) 262–266, doi: 10.1061/(ASCE)HY.1943-790 0.0 0 0 0154 .
[11] A. Defina, Two-dimensional shallow flow equations for partially dry areas, Water Resour. Res. 36 (11) (20 0 0) 3251–3264, doi:
10.1029/20 0 0WR90 0167 .
[12] C. Mügler, O. Planchon, J. Patin, S. Weill, N. Silvera, P. Richard, E. Mouche, Comparison of roughness models to simulate overland flow and tracer
transport experiments under simulated rainfall at plot scale, J. Hydrol. 402 (1-2) (2011) 25–40, doi: 10.1016/j.jhydrol.2011.02.032 .
[13] P. Costabile, C. Costanzo, F. Macchione, Comparative analysis of overland flow models using finite volume schemes, J. Hydroinf. 14 (1) (2012) 122–135,
doi:
10.2166/hydro.2011.077 .
[14] F. Simons, T. Busse, J. Hou, I. Özgen, R. Hinkelmann, A model for overland flow and associated processes within the hydroinformatics modelling system,
J. Hydroinf. 16 (2) (2014) 375–391, doi: 10.2166/hydro.2013.173 .
[15] D.P. Viero, P. Peruzzo, L. Carniello, A. Defina, Integrated mathematical modeling of hydrological and hydrodynamic response to rainfall events in rural
lowland catchments , Water Resour. Res. 50 (7) (2014) 5941–5957, doi: 10.1002/2013WR014293 .
[16] I. Özgen, K. Teuber, F. Simons, D. Liang, R. Hinkelmann, Upscaling the shallow water model with a novel roughness formulation, Environ. Earth Sci. 74
(2015) 1–16, doi:
10.1007/s12665- 015- 4726- 7 .
[17] L. Cea, M. Garrido, J. Puertas, Experimental validation of two-dimensional depth-averaged models for forecasting rainfallrunoff from precipitation data
in urban areas, J. Hydrol. 382 (1–4) (2010) 88–102, doi: 10.1016/j.jhydrol.2009.12.020 .
[18] D. Liang, I. Özgen, R. Hinkelmann, Y. Xiao, J.M. Chen, Shallow water simulation of overland flows in idealised catchments, Environ. Earth Sci. 74 (2015)
7307–7318, doi: 10.1007/s12665-015-4744-5 .
26
[19] J. Hou, Q. Liang, H. Zhang, R. Hinkelmann, An efficient unstructured MUSCL scheme for solving the 2D shallow water equations, Environ. Model. Softw.
66 (2015) 131–152, doi: 10.1016/j.envsoft.2014.12.007 .
[20] T. Dunne , W. Zhang , B.F. Aubry , Effects of rainfall, vegetation, and microtopography on infiltration and runoff, Water Resour. Res. 27 (9) (1991) 2271–
2285 .
[21] G. Blöschl , M. Sivapalan , Scale issues in hydrological modelling: a review, Hydrol. Process. 9 (1995) 251–290 .
[22] S.E. Thompson, G.G. Katul, A. Porporato, Role of microtopography in rainfall-runoff partitioning: An analysis using idealized geometry, Water Resour.
Res. 46 (7) (2010), doi: 10.1029/20 09WR0 08835 .
[23] P. Gourbesville, Data and hydroinformatics: new possibilities and challenges, J. Hydroinf. 11 (34) (2009) 330–343, doi: 10.2166/hydro.2009.143 .
[24] F. Dottori, G. Di Baldassarre, E. Todini, Detailed data is welcome, but with a pinch of salt: Accuracy, precision, and uncertainty in flood inundation
modeling, Water Resour. Res. 49 (9) (2013) 6079–6085, doi: 10.1002/wrcr.20406 .
[25] H.K. McMillan, J. Brasington, Reduced complexity strategies for modelling urban floodplain inundation, Geomorphology 90 (3-4) (2007) 226–243,
doi: 10.1016/j.geomorph.2006.10.031 .
[26] L.S. Smith, Q. Liang, Towards a generalised GPU/CPU shallow-flow modelling tool, Comput. Fluids 88 (2013) 334–343, doi: 10.1016/j.compfluid.2013.09.
018 .
[27] L. D’Alpaos, A. Defina, Mathematical modeling of tidal hydrodynamics in shallow lagoons: A review of open issues and applications to the Venice
lagoon, Comput. Geosci. 33 (4) (2007) 476–496, doi:
10.1016/j.cageo.2006.07.009 .
[28] J.-M. Hervouet , R. Samie , B. Moreau , Modelling urban areas in dam-break flood-wave numerical simulations, in: Proceedings of the International
Seminar and Workshop on Rescue Actions based on Dambreak Flood Analysis, 1, Seinâjoki, Finland, 20 0 0 .
[29] V. Guinot, S. Soares-Frazão, Flux and source term discretization in two-dimensional shallow water models with porosity on unstructured grids, Int. J.
Numer. Methods Fluids 50 (3) (2006) 309–345, doi: 10.1002/fld.1059 .
[30] S. Soares-Frazão, J. Lhomme, V. Guinot, Y. Zech, Two-dimensional shallow-water model with porosity for urban flood modelling, J. Hydraul. Res. 46 (1)
(2008) 45–64, doi:
10.1080/00221686.2008.9521842 .
[31] L. Cea, M.-E. Vázquez-Cendón, Unstructured finite volume discretization of two-dimensional depth-averaged shallow water equations with porosity,
Int. J. Numer. Methods Fluids 63 (2009) 903–930, doi: 10.1002/fld.2107 .
[32] A.S. Chen, B. Evans, S. Djordjevi
´
c, D.a. Savi
´
c, A coarse-grid approach to representing building blockage effects in 2D urban flood modelling, J. Hydrol.
426-427 (2012) 1–16, doi:
10.1016/j.jhydrol.2012.01.007 .
[33] A.S. Chen, B. Evans, S. Djordjevi
´
c, D.a. Savi
´
c, Multi-layered coarse grid modelling in 2D urban flood simulations, J. Hydrol. 470-471 (2012) 1–11, doi:
10.
1016/j.jhydrol.2012.06.022 .
[34] J.E. Schubert, B.F. Sanders, Building treatments for urban flood inundation models and implications for predictive skill and modeling efficiency, Adv.
Water Res. 41 (2012) 49–64, doi:
10.1016/j.advwatres.2012.02.012 .
[35] J. Bear , Y. Bachmat , Introduction to Modeling of Transport Phenomena in Porous Media, Kluwer Academic Publishers, Dordrecht, Boston, London, 1990 .
[36] V. Guinot, Multiple porosity shallow water models for macroscopic modelling of urban floods, Adv. Water Res. 37 (2012) 40–72, doi: 10.1016/j.
advwatres.2011.11.002 .
[37] B.F. Sanders, J.E. Schubert, H.A. Gallegos, Integral formulation of shallow-water equations with anisotropic porosity for urban flood modeling, J. Hydrol.
362 (1-2) (2008) 19–38, doi:
10.1016/j.jhydrol.2008.08.009 .
[38] M. Bruwier , P. Archambeau , S.S. Erpicum , M. Pirotton , B. Dewals , A shallow-water model with depth-dependent porosity for urban flood modeling,
in: Proceedings of the Thirty-Sixth International Association for Hydro-Environment Engineering and Research, IAHR World Congress, The Hague, the
Netherlands, 2015 .
[39] R. Hinkelmann , Efficient Numerical Methods and Information-Processing Techniques in Environment Water, Springer Verlag, Berlin, 2005 .
[40] B.R. Bird , W.E. Stewart , E.N. Lightfood , Transport Phenomena, second ed., John Wiley & Sons, New York, USA, 2007 .
[41] H.M. Nepf, Drag, turbulence, and diffusion in flow through emergent vegetation, Water Resour. Res. 35 (2) (1999) 479–489, doi: 10.1029/
1998WR90 0 069 .
[42] J. Hou, Q. Liang, F. Simons, R. Hinkelmann, A 2D well-balanced shallow flow model for unstructured grids with novel slope source term treatment,
Adv. Water Resour. 52 (2013) 107–131, doi:
10.1016/j.advwatres.2012.08.003 .
[43] C.G. Mingham , D.M. Causon , High-resolution finite-volume method for shallow water flows, J. Hydraul. Eng. 124 (June) (1998) 605–614 .
[44] Q. Liang, F. Marche, Numerical resolution of well-balanced shallow water equations with complex source terms, Adv. Water Resour. 32 (6) (2009)
873–884, doi:
10.1016/j.advwatres.2009.02.010 .
[45] W. Zhang, T.W. Cundy, Modeling of two-dimensional overland flow, Water Resour. Res. 25 (9) (1989) 2019–2035, doi: 10.1029/WR025i009p02019 .
[46] U. Razafison, S. Cordier, O. Delestre, F. Darboux, C. Lucas, F. James, A shallow water model for the numerical simulation of overland flow on surfaces
with ridges and furrows, Eur. J. Mech. B/Fluids 31 (2012) 44–52, doi:
10.1016/j.euromechflu.2011.07.002 .
[47] S. Soares-Frazão, Y. Zech, Dam-break flow through an idealised city, J. Hydraul. Res. 46 (5) (2008) 648–658, doi: 10.3826/jhr.2008.3164 .
[48] J.E. Schubert, B.F. Sanders, M.J. Smith, N.G. Wright, Unstructured mesh generation and landcover-based resistance for hydrodynamic modeling of urban
flooding, Adv. Water Resour. 31 (12) (2008) 1603–1621, doi:
10.1016/j.advwatres.2008.07.012 .
[49] B. Kim, B.F. Sanders, J.S. Famiglietti, V. Guinot, Urban flood modeling with porous shallow-water equations: A case study of model errors in the
presence of anisotropic porosity, J. Hydrol. 523 (2015) 680–692, doi: 10.1016/j.jhydrol.2015.01.059 .
[50] V. Souchere, D. King, J. Daroussin, F. Papy, a. Capillon, Effects of tillage on runoff directions: consequences on runoff contributing area within agricul-
tural catchments, J. Hydrol. 206 (3-4) (1998) 256–267, doi: 10.1016/S0022-1694(98)00103-6 .
[51] O. Planchon , M. Esteves , N. Silvera , J.-M. Lapetite , Microrelief induced by tillage: measurement and modelling of surface storage capacity, Catena 46
(2001) 141–157 .
27