DISCUSSIONS AND CLOSURES
Discussion of “Finite Volume Model for
Shallow Water Equations with Improved
Treatment of Source Terms” by Soumendra
Nath Kuiry, Kiran Pramanik, and
Dhrubajyoti Sen
February 2008, Vol. 134, No. 2, pp. 231–242.
DOI: 10.1061/ASCE0733-94292008134:2231
Alessandro Valiani
1
and Lorenzo Begnudelli
2
1
Associate Professor, Dept. of Engineering, Univ. of Ferrara, 44100 Fer-
rara, Italy. E-mail: alessandro.valiani@unife.it
2
Postgraduate Researcher, Dept. of Civil and Environmental Engineering,
Univ. of Trento, 38100 Trento, Italy.
The paper addresses a very important topic in the field of the
numerical modeling of free surface shallow flows: the treatment
of the bed slope source term in numerical methods for the solu-
tion of two-dimensional 2D shallow water equations. Using an
accurate and robust method for the treatment of this term is in fact
essential to achieving good results when dealing with complex
topography, that is, in the large majority of practical applications.
However a few points in the paper need to be discussed.
In 2D shallow water equations, the bottom slope source term
represents the horizontal component of the force acted on the
fluid by an uneven bottom. Writing the equations in differential
form, the bottom slope source term appears as S
0
= 0- ghz
b
/ x
- ghz
b
/ y
T
. In a finite volume model, the source term is inte-
grated over the computational cell, giving
S
0
d =
0-
gh
z
b
x
d -
gh
z
b
y
d
T
1
If the cell is triangular, the evaluation of the source term 1 is
trivial, as the three vertices of the cell necessarily belong to a
plane, and therefore S
0x
= z
b
/ x and S
0y
= z
b
/ y correspond to
the slopes of such plane in x and y directions, i.e., S
0x
and S
0y
are
constant over . This assumption is also made in the paper. The
source term 1 can be therefore written as
S
0
d =
0- gS
0x
hd - gS
0y
hd
T
= 0- gS
0x
V
- gS
0y
V
T
2
where V
= the volume of water stored in the cell ; using the
same notation adopted by the authors in their paper, it is simply
W
= V
. For cells that have more than 3 vertices, the evaluation
of the terms S
0x
and S
0y
is not straightforward and not unique,
since it depends on the assumptions made regarding the bottom
topography inside each cell. On this account, the technique pro-
posed by Valiani and Begnudelli 2006 is particularly useful, as
is pointed out by the authors.
From now on, we will refer to triangular cells, for simplicity
and consistency with the paper we are focusing on in this discus-
sion.
The source term is written as the product of two factors: W
and S
0x
or W
and S
0y
. Therefore, the source term does not vary
with the slope of the free surface. In fact, the volume of water W
does not vary with the free surface slope or mass conservation
would not be satisfied and bottom slopes S
0x
and S
0y
are trivially
independent from the free surface slope. Therefore, contrary to
what is stated in the paper, considering the free surface as cell-
wise constant or cell-wise linear has no effect on the evaluation of
the source term and on the final accuracy of the method.
In fact, in Eq. 22, Ah
c
is the weight of the water prism, and
the term
y
2
- y
3
h
1
+ y
3
- y
1
h
2
+ y
1
- y
2
h
3
- w
x
2A
=-
y
2
- y
3
z
1
+ y
3
- y
1
z
2
+ y
1
- y
2
z
3
2A
3
is the multiplying factor to obtain the horizontal component of
such a weight. Because h - H = z
b
, it is evident that both factors are
really independent from the free surface slope and depend on the
average free surface elevation H
c
average depth h
c
.
The cell-wise constant behavior of physical variables leads to
a final first order accuracy, while the cell-wise linear behavior
leads to a final second-order accuracy, only when we deal with
flux convective terms. It is not so when we deal with source terms:
they have a different physical and mathematical nature Hirsch
1988, so the consequences of their numerical treatment must be
properly investigated.
In fact, the method proposed by Valiani and Begnudelli
2006, where the bottom slope source term is evaluated from the
integral over the cell boundary of gh
2
/ 2, taking the free surface
elevation as a constant, is second-order accurate in space. This
has been proved both theoretically and by means of an a poste-
riori accuracy analysis Valiani and Begnudelli 2008.
Another point that is worth a short discussion is the following:
The authors write that the method of Valiani and Begnudelli
2006 could be improved by adopting the method proposed by
Komaei and Bechteler 2004. More specifically, considering the
face k, whose vertices are n
1
and n
2
, Valiani and Begnudelli
2006 evaluate the physical quantity 0.5gh
k
2
with
h
k
=
h
1
+ h
2
2
4
that is, they evaluate h
k
as the arithmetic average between the
depths h
1
and h
2
respectively at n
1
and n
2
. Instead, Komaei and
Bechteler 2004 use the following equivalent depth
h
k
=
h
1
2
+ h
1
h
2
+ h
2
2
3
5
This gives, in effect, a more accurate evaluation of the hydro-
static force on the considered face, as noted by the authors. How-
ever, it is important to note is that both methods are second-order
accurate in fact, numerical errors scale in both cases with
Ox
2
, and both methods are exactly balanced when approach-
ing still water conditions; good and nearly identical results, in-
cluding the preservation of still water conditions, can be obtained
1016 / JOURNAL OF HYDRAULIC ENGINEERING © ASCE / NOVEMBER 2009