The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1

Hyperbolic dual finite volume models for shallow water flows
in multiply-connected open channel networks○

Hidekazu Yoshioka, Doctoral Student (Research Fellow of Japan Society for the Promotion of
Science), Graduate School of Agriculture, Kyoto University, Kitashirakawa-Oiwake-cho, Sakyo-ku,
Kyoto 606-8502, Japan, yoshih@kais.kyoto-u.ac.jp
Koichi Unami, Associate Professor, Graduate School of Agriculture, Kyoto University,
Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan, unami@adm.kais.kyoto-u.ac.jp
Masayuki Fujihara, Professor, Graduate School of Agriculture, Kyoto University,
Kitashirakawa-Oiwake-cho, Sakyo-ku, Kyoto 606-8502, Japan, fujihara@kais.kyoto-u.ac.jp

Shallow water flows in multiply-connected (looped) open channel networks are encountered in various surface water
systems. Numerical simulation of such flows is one of the most important issues in environmental hydraulics and otherrelated research fields. This paper presents two numerical models of shallow water flows, dynamic wave equations andlocal inertial equations, in which junctions are implicitly and efficiently handled. Mathematical properties of the twomodels are briefly reviewed in order to understand behaviour of their solutions. Computational results of test and real
problems reveal their applicability and limitations for flows in open channel networks.
1. Introduction
Analysis of flows in surface water systems concerns a wide range ofimportant problems in environmental hydraulics and other related
research fields. Examples are flows in irrigation and drainage systems(1-2),
flash floods propagating downstream valleys(3-4) and tidal flows inestuarine systems(5-6). Macroscopic dynamics of these flows is reasonably
described with the cross-sectionally averaged one-dimensional shallow
water equations (1-D SWEs) assuming both the incompressibility ofwater and hydrostatic pressure distribution(7-8). The 1-D SWEs are acoupled system of nonlinear hyperbolic partial differential equations
having source terms. The 1-D SWEs cannot reproduce essentiallymulti-dimensional and non-hydrostatic phenomena that the 2-D and 3-Dhydrodynamic models appropriately handle(9-10), such as solitary waves,
breakings waves and oblique hydraulic jumps. Nevertheless, they have
served as one of the most effective tools for engineering purposes becauseof their efficiency. Since analytical and approximate solutions to the 1-D
SWEs are available only for limited cases, such as dam break problems inflat channels(11-12), numerical models have been used to solve them inscientific and engineering applications.
Several types of shallow water models including the 1-D SWEs havebeen presented so far. In this paper, the 1-D SWEs retaining the
momentum flux term(13) are referred to as the dynamic wave equations(DYNs) in order to distinguish them from reduced counterparts. Examplesof the reduced models are the local inertial equations (LOCs), the diffusion
wave equation (DIF) and the kinematic wave equation (KIN)(14-17). Thesemodels are derived with neglecting the temporal and/or momentum flux
terms in the momentum equation while maintaining the complete massconservation property. Although they cannot reproduce some importanttransient phenomena involving discontinuous water surface profiles thatthe 1-D SWEs deal with, their simplicity and higher efficiency are
attractive in practical applications.
Typical surface water system presents a network structure, and isregarded as a connected graph consisting of a number of junctions
connected with nodes (locally 1-D open channel networks)(18). In real
problems, shallow water flows in multiply-connected (looped) open
channel networks are frequently encountered. Development of an efficientmodeling framework for analyzing such flows is therefore an important
issue. Since a connected graph has singularities at junctions where spatial
derivatives and some hydraulic properties, such as discharge andcross-sectionally averaged velocity, are not well-defined. A key to develop
a successful numerical model is to deal with junctions efficiently andconsistently, which is achieved with the use of appropriate internalboundary conditions (IBCs) that describe balance laws of the mass andmomentum at junctions. Almost all of the existing models handle theflows in reaches and at junctions separately, resulting in the loss of
efficiency(19-21). The authors developed mathematical models of shallow
water flows in multiply-connected open channel networks that handle
junctions as implicit IBCs, which is directly incorporated into thenumerical counterparts without complicated algorithms(13, 22). Two
temporally explicit numerical models for the DYNs have been developed,
which are here refereed to be as the finite element/volume method(FEVM) that solves the continuity equation with the standard finite
element scheme and the momentum equation with a cell-centered finite
volume scheme(13), and the dual finite volume method (DFVM) thatutilizes the staggered finite volume scheme(22). The authors recently foundthat these numerical models can also be applied to the LOCs withdropping only the momentum flux term.
The purpose of this paper is to carry out theoretical and numericalanalyses to investigate applicability and limitation of the DYNs and LOCs.
These models are solved using the DFVM and applied to test and real
problems. Dam break problems and tidal flows in multiply-connected
open channel networks are considered as the test problems. As the real
problems, flows in a hydromorphic drainage system and flows in a river
network driven by flood tide are considered, both of which differ
significantly with each other in conditions. Since no direct comparison of
the DYNs and LOCs for flows in open channel networks has been carriedout, this paper contributes to further understandings of these models.
The remainder of this paper is organized as follows. Concise
introductions for the DYNs and LOCs are given in Section 2. The DFVM
is presented in Section 3. Test and real problems are carried out in
Sections 4 and 5, respectively. Section 6 provides conclusions.
2. Shallow water models
2.1 Dynamic wave equations (DYNs)
The DYNs used in this paper consist of the continuity equation¶A ¶h ¶Q
¶h ¶Q
+
-q =T
+
-q =0
¶h ¶t ¶x
¶t ¶x

1

(1)

Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1
and the momentum equation¶Q ¶ æ b Q 2 ö
æ ¶h
ö
+ ç
+ Sf ÷ = 0
÷ + gA ç¶t ¶x è A ø
¶
x
è
ø

l± =
(2)

Fb =

・Subcritical flow ( Fb2 < 1 ),・Critical flow ( Fb2 = 1 ),・Supercritical flow ( Fb2 > 1 ).

n 2Q Q

The above classification for the flow regimes is equivalent to・Subcritical flow ( F12 < b ),・Critical flow ( F12 = b ),・Supercritical flow ( F12 > b ).

(4)

A2 R 4/3

where n is the Manning’s coefficient and R is the hydraulic radius.
The DYNs comprise a system of conservative hyperbolic PDEs withsource terms, which handle both gradually and rapidly varied flows(25-26).
Momentum exchange caused by lateral inflows is not considered in (2)because only right-angled lateral inflows are assumed in this paper. The
Manning’s coefficient n and the momentum coefficient b are taken
as constant values in a domain.

If b is specified to be less than 1(28), l ± become complex undersome conditions and then the DYNs loss the hyperbolicity.
As shown in the above, the momentum coefficient b plays a central
role in determining the wave speeds. Since in practice cross-section ofchannel does not always have an idealized shape such as a rectangle, thevalue of b should be appropriately estimated. It is therefore
recommended in applications to carry out a sensitivity analysis on it.
It is also important to note influences of b on solutions to problems
with wet and dry interfaces. Among them, classical dry dam breakproblems have long been studied by many researchers in particular since
the fundamental theoretical analysis of Ritter(29) who obtained theanalytical expression of the front propagation speed of the wet and dry
interfaces(30-33). Hogg and Pritchard(33) analyzed problems with b > 1
on the basis of the characteristics theory and showed that a wet and dry
interface becomes ill-defined within a finite time if the friction slope isneglected in (2). This is due to the fact that the propagation speed of the
disturbance exceeds that of the interface at the time, which is not aphysically accepted situation.
The LOCs have the following two eigenvalues

2.2 Local inertial equations (LOCs)
The LOCs consist of the continuity equation (1) and the reduced
momentum equation
¶Q
æ ¶h
ö
+ gA ç
+ Sf ÷ = 0 ,
¶t
¶
x
è
ø

(5)

in which the momentum flux term is neglected. Being different from the
LOCs used in the literatures(15, 17), that presented in this paper deal withflows with general cross-sections.
As noted in Yoshioka et al.(27), solutions to the LOCs and DIF subject toa same boundary condition coincide if they exist. According to De
Almeida et al.(17), relative errors of water depths between the DYNs andLOCs are roughly proportional to the square of a Froude number.
Because the LOCs do not consider the momentum flux term, it is notsuited to predict flows with a large Froude number such that junctions andchannel bends are present where non-negligible horizontal exchanges ofthe momentum is important.

l± = ±

2QAgæQö± b ( b - 1) ç ÷ +
.
A
A
T
è ø

(6)

Q
±
A

Ag
.
T

(10)

2.4 Some remarks on the LOCs
This sub-section describes important remarks on the LOCs and DIF,
the latter being a doubly-nonlinear degenerate parabolic PDE(16). An
apparent difference between the LOCs and DIF is that the temporal termin the momentum equation is neglected or not. Bates and his co-workers
persist that the LOCs are superior to the DIF in the sense that an efficient
numerical model can easily be developed for the former owing to fully
retaining the temporal terms(15, 17). This is because the latter in general

Assuming a cross-sectionally uniform flow ( b = 1 ) reduces (6) to theconventional ones

l± =

Ag
,
T

which are formally derived with substituting b = 0 into (6). According
to (6) and (10), both the DYNs and LOCs are hyperbolic. As shown in(10), the wave celerity for the LOCs is slower than that of the DYNs.
Furthermore, being different from the DYNs, solutions to the LOCs do
not have critical and super-critical type flow regimes, indicating that theycannot reproduce transitional phenomena, such as hydraulic jumps andwet dam break problems.

2.3 Eigenstructures
Eigenstructures of the DYNs and LOCs are presented in thissub-section for the purpose to briefly review their differences. The DYNs
have the two eigenvalues

l± = b(9)(3)

where u is the velocity in the cross-section and V is the
cross-sectionally averaged velocity. The condition b ³ 1 holds forarbitrary shallow water flows according to the definition (3). S f is theManning’s friction slope term(24)
Sf =

F1

b -1 2
F +1
b 1

where Fb =1 in (8) is described as F1 . Local flow regimes for the
DYNs are summarized on the basis of Fb as follows:

2

1 æuö
ç ÷ dAA òA è V ø(8)which suggest the use of the “new” Froude number

where t is the time, x is the 1-D abscissa taken along the channel,
A is the cross-sectional area of flow, T is the top width of the water
surface, h is the water surface elevation, Q is the discharge, q is
the lateral inflow, g is the gravitational acceleration, b is themomentum coefficient defined as(23)

b=

ö
2
Ag æ
b -1
ç Fb =1 ±( Fb =1 ) + 1 ÷÷ ,
T çè
b
ø

(7)

The eigenvalues l ± in (6) are rewritten using the “conventional”
Froude number Fb =1 as
2Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1
straightforward because it is achieved simply neglecting the momentum
flux term of (2). Spatial discretization for the other terms is carried out in asimilar manner. Since the LOCs do not have the momentum flux term in
its formulation, no momentum flux evaluation scheme is necessary inorder to solve them in locally 1-D open channel networks.
Temporal integration of the DFVM is carried out using the classicalfourth-order Runge-Kutta method. Since this temporal integration method
is of an explicit type, the time increment Dt has to be chosen as asufficiently small value(44) such that

requires the use of a very small time increment for temporally explicit
method and a time-consuming iterative algorithm for implicit method.
Another remark to be noted is a relationship between the LOCs and
DIF in an asymptotic limit. According to the theory of Chapman-Enskog
expansion for singular hyperbolic PDEs(34), formally introducing a small
parameter 0 < e = 1 in (5) yields

e

¶Q
æ ¶h
ö
+ gA ç
+ Sf ÷ = 0 .
¶t
¶
x
è
ø

(11)

Taking the limit of e ® +0 in (11) leads to the DIF, a parabolic system
of PDEs consisting of the continuity equation (1) and¶h+ Sf = 0 .
¶x

Dx
Ag
> b V + b ( b - 1)V 2 +
,
Dt
T

(13)

where Dx is the elemental length is at least satisfied in an entire
computational domain.

(12)

Application of an asymptotic-preserving scheme(35) to (1) and (11) with asufficiently small e yields the approximate solution of the DIF withnumerically retaining the temporal term of (11), which is far efficientcompared with directly solving the DIF. This research topic is beyond thescope of this paper and will be addressed in future researches.

4. Test problems
Since the DFVM has been extensively verified with test problems inopen channels, this paper focuses mainly on those of in open channel
networks. An exception is the dam break problems firstly carried out inthis section to see influences of b on solutions to the DYNs.

3. Dual finite volume method (DFVM)
In this paper, the DYNs and LOCs are solved using the DFVM
(Appendix), a recently developed simple and versatile numerical method
for shallow water flows both in single open channels and in open channel
networks that have channel bends and multiply-connected structures. The
DFVM applies a node-centered finite volume scheme to the continuity
equation (1) and an upwind, cell-centered finite volume scheme to themomentum equation (2), respectively. The water surface elevation hand the discharge Q serve as the dependent variables, which are
distributed to the dual cells and regular cells respectively. The dual cell isthe 1-D counterpart of the Voronoi diagram for multiple dimensions.
Upwind algorithms applied to the momentum flux and source terms
stabilize numerical solutions without explicitly adding artificial
dissipation terms. The DFVM does not require the use of water surface
reconstruction algorithms to preserve a still water condition because it
utilizes the water surface elevation h as a dependent variable. Local
mass conservation at junctions is achieved as the implicit IBCs that do not
rely on complicated numerical algorithms.
In this paper, the DFVM is slightly modified so that the upwinding of
the momentum flux term in (2) is carried out with the new Froudenumber Fb not with the conventional one F1 . In the original DFVM,
local flow regimes are identified from the value of F1 , which is notconsistent with the eigenstructures of the DYNs as indicated in theprevious section. According to preliminary numerical tests not presented
in this paper, the modification does not significantly improve solution
profiles; however, it certainly enhances theoretical foundations of theDFVM. In addition to the above modification, a new momentum flux
evaluation scheme to more accurately simulate flows around junctions is
implemented in the spatial discretization procedure(36). Incorporation of
the scheme into the DFVM significantly improves the accuracy to predictthe discharge ratios for the flows downstream of a junction. The accuracy
in some cases is comparable to that of a horizontally 2-D shallow watermodel, demonstrating its validity(37-38). Performances of the DFVM havebeen verified with a series of benchmark tests of subcritical, supercriticaland transcritical flows(39-42), clearly showing that it is as accurate as theother recently developed numerical models. The FEVM, a finite/element
counter part of the DFVM, has also successfully been applied to a wide
range of problems(13, 35, 43).
Application of the DFVM to the LOCs presented in this paper is

4.1 Dam break problems in a single channel
Dam break problems in an open channel serve as fundamental test
cases to examine ability of a numerical model for the DYNs to resolve
upstream and downstream wave propagations involving transitions andshocks. Here dam break problems in a straight, flat channel areconsidered. A 2,000 (m) length open channel is considered as the domain
W = ( 0,2000 ) (m). Cross-section of the channel is a rectangle with thewidth of 1 (m). The initial water elevation h (m) is specified as
ïì hU
ïî hD

h =í

( x £ 1000 )( x > 1000)(14)where hU (m) and hD (m) denote the initial upstream anddownstream water depths, respectively. The initial discharge is Q = 0(m3/s) in the entire W . The initial upstream water depth hU is fixed to
10 (m), while the initial downstream water depths hD are considered
three cases: 5.0 (m) (Case DBA), 0.1 (m) (Case DBB) and 0.0 (m) (Cases
DBC and DBD). Cases DBA and DBB are the problems in a frictionlesschannel, which intend to reveal dependence of the water surface profiles
on b . On the other hand, in Cases DBB and DBC, non-zero Manning’s
coefficient n are specified as 0.005 (s/m1/3) in Case DBC and 0.020(s/m1/3) in Case DBD, respectively, both of which are realistic values inexperimental and fields situations. The upstream and downstream
boundaries of W are solid walls where flows are reflected. W is
uniformly discretized into 2,000 regular cells so that numerical solutions
sufficiently converge. The time increment Dt is fixed to 0.0015 (s).
The computed water surface profiles are illustrated in Figs. 1 through 4
for each case with different values of b , namely, b = 1.00 , 1.05, 1.10,
1.15 and 1.20. The label “DBA-1.00” represents the water surface profilein Case DBA for the value of b = 1.00 , and so is for the others. Figs. 1
and 4 show that the propagation speeds of the shocks and wet and dry
interfaces increase as b increases. Figs. 1 and 2 show that water depth
immediately upstream of the shocks decrease smaller as b increases.
Figs. 3 and 4 indicate that the water surface profiles are apparently lesssensitive to b in Cases DBC and DBD than those in the wet channel
bed counterparts, but its influence is more significant for the smaller n .
Although not presented here, specifying n = 0 (s/m1/3) in a dry dam
break problem results in a numerical failure being consistent with thetheoretical results of Hogg and Pritchard(33).
3

Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1

Fig. 1 Computed water surface profiles for Case DBA ((a) t = 15 (s), (b) t = 30 (s) and (c) t = 45 (s))

Fig. 2 Computed water surface profiles for Case DBB ((a) t = 15 (s), (b) t = 30 (s) and (c) t = 45 (s))

Fig. 3 Computed water surface profiles for Case DBC ((a) t = 15 (s), (b) t = 30 (s) and (c) t = 45 (s))

Fig. 4 Computed water surface profiles for Case DBD ((a) t = 15 (s), (b) t = 30 (s) and (c) t = 45 (s))
4.2 Dry dam break problem in an open channel network
A dry dam break problem in a hypothetical multiply-connected open
channel network (Fig. 5) is carried out with the DYNs. This test problem
is challenging because it requires a consistent junction treatment techniqueas well as a stable discretization method for rapidly varying transcriticalflows involving wet and dry interfaces, non-rectangular cross-sections andchannel bends. Only a few attempts at cross-sectionally averaged
modeling have been made for the problems in open channel networks.
4

Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1

Fig. 5 Open channel network for test problems

Fig. 7 Computational domain for the first real problem

Fig. 6 Water surface profiles in the open channel network for the DYNs
presented at the time interval of 5 (s)The reaches in the channel network are equal in terms of length and bed
slope, which are 10 (m) and 0.01, respectively. Cross-section of thechannel is a rectangle with the width of 10 (m) in the reach A-H and atriangle with the side slope of 1:1 in the others. Manning’s coefficient n
is set as 0.040 (s/m1/3) for the reaches B-C-D and 0.030 (s/m1/3) for theothers. The momentum coefficient b is set as 1.1. A vertical wall
serving as a dam is installed at the middle of reach A-B, namely H,
dividing the channel network into an upstream reservoir (A-H) and adownstream multiply-connected channel network. The initial water depth
immediately upstream of the dam is set as 0.6 (m) and the downstream
channel network is initially dry. Both the upstream and downstream-ends
are solid walls. The dam is instantaneously removed at the initial time
t = 0 (s). Each reach is divided into 200 regular cells. The timeincrement Dt is set as 0.001 (s).
Fig. 6 shows the computed water surface profiles in the domain. A
surge from the reservoir separates at B. The divided surges converge at D
and hit the downstream wall at E, generating a receding bore propagating
upstream. There are also other receding bores starting at bending points C
and F, as observed in some experiments(45). Increasing spatial andtemporal resolutions do not significantly alter the computational results,
showing validity of the numerical simulation carried out here.
The LOCs has also been applied to this problem; however, thecomputed water surface profile diverged, highly oscillating around thedam H. This is considered due to the instability of the DFVM for theLOCs, indicating necessity of improvements.

Fig. 8 Bed elevation distribution in the domain

Fig. 9 A sketch of cross-section for the first real problem
flows in a hypothetical flat, multiply-connected open channel network.
The domain has a same network structure with the locally 1-D open
channel network used in the previous test problem. Length of each reachis 1,000 (m). Cross-section of the channel is a rectangle with the width of
25 (m) for all the reaches. The momentum coefficient b is set as 1.1.
The boundary E is a solid wall. The boundary A is treated as an open
boundary at which the water surface elevation h A is specified directly
using a sinusoid function asæ 2p t ö÷.
è 43,200 ø

h A = 5 + sin ç

4.2 Tidal flows in an open channel network
The DYNs and LOCs are next applied to numerical simulation of tidal
5(15)

Initial condition is set as a still water with the water depth of 5 (m). Each
reach is divided into 200 regular cells. The time increment Dt is set as
0.5 (s). Terminal time of the computation is 43,200 (s).
Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1
Maximum relative errors between the computed water depths for the
DYNs and LOCs at the key nodes B, C, D, E and F are 9.2×10-4, 1.1×10-3,
1.2×10-3, 1.2×10-3 and 1.2×10-3, respectively, demonstrating that theresults of the two models are almost identical. This is because the Froudenumber of the flow is sufficiently small such that the use of the LOCs isjustified. The results indicate applicability of the LOCs to numericalsimulation of tidal flows in multiply-connected open channel networks
with relatively wide channel widths.
5. Real problems
For the purpose of making further comparisons of the DYNs and LOCs,
numerical simulations of two significantly different real shallow waterflows in multiply-connected open channel networks are carried out.
5.1 Flows in a hydromorphic drainage system
Applicability of the DYNs and LOCs to real problems is firstly
assessed with surface water flows in a multiply-connected open channel
network draining hydromorphic farmlands in the Guinea savanna
agro-ecological zone of Ghana(13, 46). The main water input into thechannel network is direct rainfall and lateral groundwater seepage fromthe surrounding rice fields. Fig. 7 illustrates a plane view of the open
channel network as the computational domain with the mesh and key
nodes alphabetically labeled from A through P. The open channel network
consists of a main loop (A-B-C-D-E-F-A) running in the rice fields, adownstream gully (K-L-M-N-O-J-P) and steep cliffs connecting the main
loop and the gully (F-L, A-M, G-N, H-O and I-J). The node D is thehighest point in the channel network, K and P the boundaries, A, F, G, H, J,
L, M, N and O are junctions. The node P is the downstream-end of thedomain where a free-overflow natural weir is installed. Fig. 8 shows bed
elevation distribution in the domain where that of the downstream-end P
is set as the normal elevation 0 (m). Bed slopes in the domain are very
irregular with the maximum value of 1/2.14 in a gully. Cross-section of
the channel is non-rectangular and non-prismatic, which is governed by
the five parameters H , Wl , Wr , Sl and S r as shown in Fig. 9(13).
Total numbers of the regular cells and the nodes in the computational
domain are 228 and 224, respectively. Yoshioka et al.(18, 47) computed
mean residence time distribution of a conservative solute injected into thedomain on the basis of Kolmogorov’s backward equations, estimating the
maximum mean residence time of the solute in the domain as 2,200 (s) at
a no-rainfall steady state and as 500 (s) that of during the rainfall event
focused on later, respectively.
The momentum coefficient b and the Manning’s coefficient n in
the domain have been estimated as 1.10 and 0.015 (s/m1/3), respectively(13).
Surface water flows in the channel network exhibit highly complicatednatures including both subcritical and supercritical flows due to its
irregular topography. The time increment is set as Dt = 0.05 (s). Other
computational conditions are same with those in Unami and Alam(13).
Fig. 10 presents computational results of the outflow hydrographs at thedownstream-end P for the DYNs and LOCs during a rainfall event in arainy season that caused severe ﬂoods in downstream area of the
farmlands. The outflow discharge of the LOCs is gently varying
compared with that of the DYNs. Maximum relative error between the
discharges for the models is more than 0.4.
Another set of numerical simulations for the DYNs are also carried out
under the same initial and boundary conditions with varying the values ofthe momentum coefficient b for the purpose of a brief sensitivity
analysis. Computational results of the relative errors for the outflow
discharge at the downstream-end P are shown in Fig.11. The label“Relative error-1.00” represents the relative error of the outflow

Fig. 10 Computed outflow discharges at P for the DYNs and LOCs
and the relative errors between themFig. 11 Computed outflow discharges at P for the DYNs withdifferent values of the momentum coefficient bdischarges between that of the estimated value b = 1.10 and that ofb = 1.00 , and so is for the others. Maximum value of the relative error isless than 0.08 for these computational cases, showing that varying b
within the presented range does not significantly alter the hydrograph.
5.2 Flows in a river network driven by a flood tide
The DYNs and LOCs are next applied to numerical simulation of flows
driven by a flood tide in a tidal river network surrounding a polder near
the coast of Bay of Bengal, Bangladesh(48). Fig. 12 illustrates a plane viewof the river network as the computational domain with labeled key nodes.
The river network is connected with Bay of Bengal at B-0. The boundary
B-1 is connected with one of the downstream tributaries of the MeghnaRiver. A polder surrounded by tidal river is vulnerable to flood tides that
may cause a severe flooding events associating salt damages(49-50). In
these coastal areas, such a serious event is caused by the arrival of atropical typhoon. Shallow water models can serve as the foundation ofeffective analysis of this kind of disasters.
The river network is discretized into a mesh with 398 regular cells and
394 nodes. The normal water level is set as 0 (m) above the sea level.
Channel bed elevation in the large river network surrounding the polder is
set as -10 (m), while those for the small channels in the polder are
distributed between -4 (m) and 0 (m). A still water with the surface water
elevation of 0 (m) in the entire domain is set as the initial condition. The
two boundaries B-0 and B-1 are open boundaries at which the external
6

Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1
water surface elevation h E for a hypothetical flood tide event is
specified as(48)æt86,400 - t ö ö,÷ ÷ (m).
43,
200
43,200 ø øèæh E = max ç 0,5min çè(16)

The discharge at the open boundaries B-0 and B-1 are specified on thebasis of h E as
QE = TB-0 or B-1 hB-0 or B-1 - h E

2 g h B-0 or B-1 - hE

(m3/s).

(17)

The other boundaries are treated as solid walls. The time increment is set
as Dt = 1 (s).
Figs. 13 and 14 show the computed water surface elevations in the river
network at the peak ( t = 43,200 (s)) and at the end ( t = 86,400 (s)) of
the flood tide for the DYNs and LOCs, respectively. Relative errors
between the water surface elevations at the peak and end of the flood tideare 1.1×10-2 and 7.7×10-3, respectively, showing that the LOCs givecomparative computational results with the DYNs at the peak. The
computed water surface elevations for the LOCs are less than those of the
DYNs, indicating that application of the former to analysis of a flood tide
event results in an underestimation of its risk.
Although the simulation has successfully been carried out without
numerical failures, it should be emphasized that the DFVM presented inthis paper computes flows only in the channels. Overtopping flows fromthe channels to the polders are therefore not taken into account. Morerealistic, detailed numerical simulations of flood tides require the use of a
sophisticated coupled modelling framework of water flows both incross-sectionally 1-D open channel flows and horizontally 2-D overland
flows(51-52). An FEVM for solving the 2-D SWEs currently underdevelopment can be used for this purpose(53).

Fig. 12 Computational domain for tidal flows in an estuary

6. Conclusions
Two shallow water models, the DYNs and LOCs, are presented andapplied to test and real problems in order to investigate their applicability
and limitations. Several mathematical properties of the DYNs and LOCs
were briefly reviewed mainly focusing on their eigenstructures.
Computational results of the test and real problems for the DYNs with the
DFVM demonstrated their versatility, well reproducing the numericalsolutions without numerical failures under the severe computational flows.
The computational results also suggested high applicability of the LOCs
to analysis of the flows in multiply-connected tidally-driven river
networks where the Froude number is sufficiently small. The results werecomparable to those of with the DYNs. On the other hand, thecomputational results of the dam break problems showed clearly that thepresent LOCs should not be applied to simulate these kinds of flows. In
conclusions, the LOCs give reasonable results for the flows with small
Froude numbers where the momentum exchanges are negligible.
As demonstrated in this paper, the DFVM, a fully explicit numericalmodel for the DYNs and LOCs, has successfully been applied tonumerical simulations under various flow conditions. Computational
performance of the DFVM can be improved if a semi-implicit algorithm
is incorporated into it in order to remove the restriction for the timeincrement. Then the improvements should be made so that its algorithmic
simplicity is not degraded. There remain many important research topics
to be addressed in order to develop a more promising model for analyzing
a wide range of shallow water flows, such as coupling of the DYNs withmulti-dimensional shallow water models or with groundwater flow
models(54-55). The DFVM used in this paper will be applied to risk analysis
of a catastrophic dam break flash flood in Japan caused by a huge
earthquake, and the analysis results will be presented in elsewhere.

Fig. 13 Computed water surface profiles for the DYNs and LOCs at
the peak of the flood tide

Fig. 14 Computed water surface profiles for the DYNs and LOCs at
the end of the flood tide

7

Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1
where å q represents the discharge contributed from the lateral
inflows and the boundary conditions.

Acknowledgement
This paper is financially supported by the JSPS under grant 25･
2731.
Appendix
Spatial discretization algorithm for the DFVM used in this paper is
described in this appendix.

A.3 Momentum equation
The upwind FVM(13) based on the local Froude number is applied tothe momentum equation (2). The cell-centered finite volume formulationof the momentum equation in W k leads to

A.1 Computational mesh
A couple of staggered computational meshes, a regular mesh and a dual
mesh, are used in the DFVM. The continuity equation (1) is solved on the
dual mesh, while the momentum equation (2) on the regular mesh, so as
to avoid conflicting numbers of equations and unknowns. The domain
W is first divided into a regular mesh consisting of regular cells boundedby two nodes, so that a junction falls on one of the nodes. The regularcells and the nodes are indexed with the natural numbers. The total
numbers of regular cells and nodes are denoted by Nc and N n ,respectively. The kth regular cell is denoted by W k . The length of W k
is represented by lk . The two nodes bounding both sides of W k aredenoted by the j ( k ,1) th node and the j ( k ,2 ) th node. The abscissain W k is directed to the j ( k ,2 ) th node. The number of regular cells
connected to a generic ith node is denoted by v ( i ) . The jth regular cellconnected to the ith node is referred to as the k ( i, j ) th regular cell
Wk (i , j ) . There are two nodes that bound Wk (i , j ) . One is the ith node,and the other is referred to as the m ( i, j ) th node. In Wk (i , j ) , thedirection of the abscissa is identified with the parameter s i , j , whichequals 1 when x is directed to the m ( i, j ) th node, and otherwise to -1.
A dual mesh is generated from the regular mesh. Following themulti-dimensional analogue of Mishev(56), define the ith dual cell Si as{}

S i = x xi - x < xm ( i, j ) - x for 1 £ j £ n (i)

lk

Vk =

Si

i, j

i, j

Si

x - xis i , jlk (i , j )(Ti , j ,1

- Ti , j ,0 )

æ

(18)

wk = max ç 1 è

ö
1
,0 ÷
2
Frk
ø

(26)

where Frk2 is the square of the cell local Froude number Fb in (9)computed using b , g , Tk ,DS , Qk and Ak , DS , and the subscripts
US and DS indicate values at the upstream node and the downstream
node, respectively, in W k . The velocity Vk is taken to be 0 when
Ak ,DS < e . Finally, the cell flux Fk is determined as
Fk = b k QkVk .

(27)

The flux on the cell interface ¶W k is evaluated considering connections
between cells and the momentum balance principles(13, 35). The set ofindices of the regular cells, the downstream (upstream) nodes of whichfall on the upstream (downstream) node of W k , is denoted as U k( Dk ). The flux at the upstream cell interface ¶Wk ,US Ì ¶Wk isprescribed as

(19)

F ¶W

k ,US

= å rk max ( cos q k ,k ,0 ) Fk

(28)

k ÎU kwhere rk is the discharge ratio defined analogues to that of Yoshioka etal.(36) as

(20)

æ
ö
rk = Qk çç å Qk ÷÷è k ÎDkø

where Ti , j ,0 and Ti , j ,1 are the values of T at the ith node and at them ( i, j ) th node in Wk (i , j ) , respectively. The first term on the left-handside of (19) is evaluated as

-1

(29)

(21)

which equals to 0 if its denominator vanishes and q k ,k is the angle
between c k and ck . The discharge ratio rk satisfies the partition ofunity property(30)
å rk = 1 .

æ 1 n (i )öçç å lk (i , j ) ( 3Ti , j ,0 + Ti , j ,1 ) ÷÷ (22)
è 8 j =1
øThe flux at the downstream cell interface ¶Wk ,DS Ì ¶Wk is identical tothe cell flux

æ 1 n (i )
ö dh
¶h
T
d
x
=
çç å lk (i , j ) ( 3Ti , j ,0 + Ti , j ,1 ) ÷÷ i ,òSi ¶tè 8 j =1
ø dt

leading to
n i
ö
dh i æ ( )= ç - ås i, j Qk (i , j ) + å q ÷÷dt çè j =1
ø(25)

and the weight

Assuming a linear variation of the top width T in Wk (i , j ) yields thefollowing relationship
T = Ti , j ,0 +(24)

Ak = (1 - wk ) Ak ,DS + wk Ak ,USn (i )

j =1

Qk
Ak

with the cell cross-sectional area

A.2 Continuity equation
The continuity equation (1) is discretized on a dual mesh. The
cell-vertex finite volume formulation of the continuity equation in the
dual cell Si is
¶h

(23)

where W k is the interface of ¶W k .
Flux evaluation in the momentum equation (2) is carried out with asimple upwind discretization. For each generic regular cell W k , the node,of the two nodes bounding W k , to which the flow is directed is referred
to as the downstream node, and the other is referred to as the upstream
node. The vector starting from the upstream node and ending at thedownstream node is denoted by c k . The cell flux Fk for the regularcell W k is determined using the local Froude number as a weight.
When the downstream node in W k is wet, Ak ,DS ³ e for a small
threshold value e , and the cell cross-sectional area Ak and the cellcross-sectionally averaged velocity Vk are calculated aswhere xi and xm (i , j ) represent x at the ith node and at them ( i, j ) th node, respectively. The dual mesh consists of N n dual cells.
The cell interface between Si and S m (i , j ) is denoted by G i, j . The
water elevation h is attributed to the dual cells, and the discharge Q
to the regular cells. The discretized h in Si is denoted by hi , and thediscretized Q in W k by Qk .

ò T ¶t dx + å s Qk ( ) = ò qdx .

dQk
æ ¶h
ö
+ [ F ]¶W = ò gA ç - S f ÷ dx
k
Wk
dt
¶
x
è
ø

k ÎDk

8

Copyright © 2013 by JSFM

= Fk

(31)

= Fk - min ( cosq k ,k ,0 ) Fk

(32)

F ¶W

k ,DS

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1(4) Yoshioka, H., Unami, K. and Fujihara, M., “A dual grid based finite
volume model for dam break flows in multiply connected open
channel networks with general cross-sections”, Proc. the 62nd
Japan National Congress for Theoretical and Applied Mechanics,
2013a, OS19-09, 2pp.
(5) Bushman, F.A., Hoitink, A.J.F., Van Der Vegt, M. and Hoekstra, P.,“Subtidal flow division at a shallow tidal junction”, Water Resour.
Res., 46 (2010), W12515, 12pp.
(6) Sassi, M.G., Hoitink, A.J.F., De Brye, Vermeulen, B. andDeleersnijder, E., “Tidal impact on the division of river dischargeover distributary channels in the Mahakam Delta”, Ocean Dyn.,
61(12) (2011), pp. 2211-2228.
(7) Cunge, J.A., Holly Jr., F.M. and Verwey, A., Practical Aspects ofComputational River Hydraulics, 1980, Pitman, pp. 7-52.
(8) Szymkiewicz, R., Numerical Modeling in Open Channel
Hydraulics, Springer, 2010, pp.263-300.
(9) Yang, C., Lin, B., Jiang, C. and Liu, Y., “Predicting near-fielddam-break flow and impact force using a 3D model”, J. Hydraul.
Res., 48(6) (2011), pp. 784-792.
(10) Shinha, S., Liu, X. and Garcia, M.H., “Three-dimensionalhydrodynamic modeling of the Chicago River, Illinois”, Environ.
Fluid Mech., 12(5) 2012, pp. 471-494.
(11) Stoker, J.J., Water waves. Interscience publishing, New York, 1957,
pp.291-341.
(12) Delestre, O., Lucas, C., Ksinant, P.A., Darboux, F., Laguerre, C., Vo,
T.N.T., James, F. and Cordier, S., “SWASHES: a compilation ofshallow water analytic solutions for hydraulic and environmental
studies”, Int. J. Numer. Methods Fluids, 72(3) (2013), pp.269-320.
(13) Unami, K. and Alam, A.H.M.B., “Concurrent use of finite element
and finite volume methods for shallow water equations in locally
1-D channel networks”, Int. J. Numer. Methods Fluids, 69(2) (2012),
pp. 255-272.
(14) Yen, B.C. and Tsai, C.W.S., “On noninertia versus diffusion wave
in flood routing”, J. Hydrol., 244(1-2) (2001), pp. 97-104.
(15) Bates, P.D., Horritt, M.S. and Fewtrell, T.J., “A simple inertial
formulation of the shallow water equations for efficient
two-dimensional flood inundation modelling”, J. Hydrol., 387(1-2)(2010), pp. 33-45.
(16) Santillana, M. and Dawson, C., “A numerical approach to study theproperties of solutions of the diffusive wave approximation of theshallow water equations”, Comput. Geosci., 14(1) (2010), pp.
31-53.
(17) De Almeida, G.A.M. and Bates, P., “Applicability of the local
inertial approximation of the shallow water equations to floodmodeling”, Water Resour. Res., 49 (2013), 12pp.
(18) Yoshioka, H., Unami, K. and Kawachi, T., “Partial differential
equation model for spatially distributed statistics of contaminantparticles in locally one-dimensional open channel networks”, Proc.
Tenth International Conference on Hydroscience and Engineering,(2012), 24pp.
(19) Colombo, R.M., Guerra, G., Herty, M., and Schleper, V. “Optimal
control in networks of pipes and canals”, SIAM, J. Control Optim.,
48(3) (2009), pp. 2032-2050.
(20) Marigo, A., “Entropic solutions for irrigation networks”, SIAM, J.
Appl. Math., 70(5) (2010), pp. 1711-1735.
(21) Herty, M. and Seaid, M., “Assessment of coupling conditions inwater way intersections”, Int. J. Numer. Methods Fluids, 71(11)(2012), pp. 1438-1460.
(22) Yoshioka, H. and Unami, K., “A dual grid based finite volume

with the exception that
F ¶W

k ,DS

when the downstream node is connected to exactly two regular cells.
Finally, the second term of the left-hand side of (23) is evaluated using
F ¶W
and F ¶W
with considering the directions of the flows as
k ,US
k ,DS
well as that of the x abscissa. The momentum flux evaluation scheme
presented is physically relevant in the sense that it considers the contact
angles of the channels meeting at a junction and guarantees non-increaseof the momentum variation(57). Application of an improper, momentum
variation-increasing scheme to the DFVM leads to unphysical surface
water profiles(36, 58).
The source terms in the right-hand side of (23) are discretized using an
upwind algorithm. The first term of the right-hand side of (23) isevaluated asæ ¶h ö

ò gA çè - ¶x ÷ø dx = - gl A (hj ( ) - hj ( ) )(33)

ìAk , US( Frk 2 ³ 1)
ï
Ak = í A.
+ Ak ,j ( k ,2)
k ,j ( k ,1)
ï
Frk 2 <1)
(
2
î

(34)

k

Wk

k

k ,2

k ,1

with

The second term of the right-hand side of (23), the friction slope term, isdiscretized asò gA ( - S ) dx = - gl A Skknk2Qk Qk P k

4/ 3

f

Wk

(35)

f ,k

with
Sf , k =

Ak

(36)

10/3

and
ì
Pk , US( Frk 2 ³ 1)
ï
Pk = í P.
+ Pk ,j ( k ,2)
k ,j ( k ,1)
ï
Frk 2 <1)
(
2
î

(37)

However, S f is taken as 0 when Ak < e . Finally, each
explicitly computed as

(

(

dQk 1
=
- [ F ]¶W - g Ak hj ( k , 2) - hj ( k ,1) + lk S f , k
k
dt
lk

)) .

dQk
is
dt

(38)

The system of ordinary differential equations (22) and (38) is solved
using the Runge-Kutta method with appropriately specified initial andboundary conditions.
Bibliography(1) Trancoso, A.R., Braunschweig, F., Leitao, P.C., Obermann, M. andNeves, R., “An advanced modeling tool for simulating complex
river systems”, Sci. Total Environ., 407(8) (2009). pp. 304-3016.
(2) Van Thang, P., Chopard, B., Lefèvre, L., Ondo, D.A. and Mendes,E., “Study of the 1D lattice Boltzmann shallow water equation andits coupling to build a canal network”, J. Comput. Phys., 229(19)(2010), pp. 7373-7400.
(3) Carrivick, J.L., Manville, V., Graettinger, A. and Cronin, S.J.,“Coupled fluid dynamics-sediment transport modelling of a Crater
Lake break-out lahar: Mt. Ruapehu, New Zealand”, J. Hydrol.,388(3-4) (2010), pp. 399-413.
9

Copyright © 2013 by JSFM

The 27th Computational Fluid Dynamics Symposium
Paper No. B07-1(42) Aricò, C., Sinagra, M. and Tucciarelli, T., “Anisotropic potential of
velocity ﬁelds in real ﬂuids: Application to the MAST solution ofshallow water equations”, Adv. Water Resour., 62(Part A) (2013),
pp.13-26.
(43) Kinjo, N., Yoshioka, H., Unami, K. and Fujihara, M., “Simulationof water flows in an agricultural drainage system with hydraulicstructures”, Proc. 2013 IAHR World Congress, (2013), 12pp.
(44) Toro, E.F., Riemann Solvers and Numerical Methods for Fluid
Dynamics, Springer, 2010, pp. 163-212.
(45) Fazãro, S.S. and Zech, Y., “Effects of a sharp bend on dam-breakflow”, Proc. the 28th IAHR congress, (1999), 7pp.
(46) Unami, K., Kawachi, T., Kranjac-Berisavljevic, G.., Abagale, F.K.,
Maeda, S. and Takeuchi, J., “Case study: Hydraulic modeling of
runoff processes in Ghanaian inland valleys”, J. Hydraul. Eng.,
135(7) (2009), pp. 539-553.
(47) Yoshioka, H., Unami, K. and Kawachi, T., “An upwind finite
element scheme for stochastic contaminant particle dynamics inlocally one-dimensional open channel networks”, Proc. PAWEES
2011 International Student’s Conference, (2011), 16 pp.
(48) Ran, E., Yoshioka, H., Unami, K., and Fujihara, M.,“Computational fluid mechanics-based investigation of operational
strategies for drainage gates installed in a polder”, Proc. the Annual
Meeting of Applied Hydraulics, (2012), pp. 104-105. (in Japanese)(49) Debsarma, S.K., “Simulations of storm surges in the Bay of
Bengal”, Mar. Geod., 32(2) (2009), pp. 178-198.
(50) Lewis, M., Bates, P., Horsburgh, K., Neal, J. and Schumann, G., “A
storm surge inundation model of the northern Bay of Bengal using
publicly available data”, Q. J. R. Meteorol. Soc. 139(671) (2013), pp.
358-369.
(51) Chen, Y., Wang, Z., Liu, Z. and Zhu, D. (2012): “1D–2D coupled
numerical model for shallow-water flows”, J. Hydraul. Eng., 138(2)(2012), pp. 122–132.
(52) Morales-Hernández, M., García-Navarroa, P., Burguete, J. andBrufau, P., “A conservative strategy to couple 1D and 2D models
for shallow water flow simulation”, Comput. Fluids, 81(20) (2013),
pp. 26-44.
(53) Yoshioka, H., Takagi, K., Unami, K. and Fujihara, M., “Ahorizontally 2-D shallow water model utilizing both the finite
element and the finite volume techniques”, Proc. the 70th JapaneseSociety of Irrigation, Drainage and Rural Engineering, Kyoto
Division, (2013), 2pp. (in Japanese, to appear)(54) He, Z., Wu, W. and Wang, S.S.Y., “Coupled finite volume model for
2D surface and 3D subsurface flows”, J. Hydraul. Eng., 13(9)(2007), pp. 835-845.
(55) Xin, P., Yuan, L.R., Li, L. and Barry, D.A., “Tidally driven
multiscale pore water flow in a creek-marsh system”, Water Resour.
Res., 47(W07534) (2011), 19pp.
(56) Mishev, I.D., “Finite volume methods on Voronoi meshes”, Numer.
Methods Partial Differ. Equat., 14(2) (1998), pp. 193-212.
(57) Yoshioka, H., Unami, K. and Fujihara, M., “Relevant internalboundary conditions for the partial differential equations In locally
one-dimensional open channel networks”, Proc. the Annual
Meeting of Applied Hydraulics in 2013, (2013), 6pp. (to appear)(58) Ishida, K., Yangyuoru, M., Unami, K. and Kawachi, T. “Applicationof shallow water equations to analyze runoff processes in hilly
farmlands”, Paddy Water Environ., 9(4) (2011), pp. 393-401.

scheme for analysis of shallow water flows in locally
one-dimensional open channel network flows”, Proc. the 69th
JSDIRE, Kyoto Division, (2012), pp. 114-115. (in Japanese)(23) Costabile, P. and Macchione, F., “Analysis of one-dimensionalmodelling for flood routing in compound channels”, Water Resour.
Manage., 26(5) 2012, pp. 1065-1087.
(24) Yen, B.C., “Open channel flow resistance”, J. Hydraul. Eng, 128(1)(2002), pp. 20-39.
(25) Lai, C., Baltzer, R.A. and Schaffranek, R.W., “Conservation-form
equations of unsteady open-channel flow”, J. Hydraul. Res., 40(5)(2002), pp. 567-578.
(26) Gąsiorowski, D., “Balance errors generated by numerical diffusionin the solution of non-linear open channel flow equations”, J.
Hydrol., 476 (2013), pp. 384-394.
(27) Yoshioka, H., Kinjo, N., Wakazono, A., Unami, K. and Fujihara, M.,“Application of two shallow water models to steady flow analysis ina vegetated agricultural drainage canal”, J. Rainwater Catch. Syst.
(under review)(28) Abedini, M.J. and Hashemi, M.R., “Effect of convective term
suppression in numerical simulation of trans-critical open channel
flow”, Iranian J. Sci. Tech. Trans. B. Engrg., 31(B1), pp. 85-96.
(29) Ritter, A., “Die Fortpflanzung der Wasserwellen”, Z. Vereines
Deutsch. Ing., 36(33) (1982), pp. 947-954. (in German)(30) Whitham, G.B., “The effects of hydraulic resistance in the
dam-break problem”, Proc. R. Soc. Lond. A., 227 (1955), pp.
399-407(31) Chanson, H., “Analytical Solutions of Laminar and Turbulent Dam
Break Wave”, Proc. of River Flow, 2006, pp. 465-474.
(32) Chen, Y.C., Wu, B. and Wang, B., “Similarity solution of dam-breakflow on horizontal frictionless channel”, J. Hydraul. Res., 49(3)(2012), pp. 384-387.
(33) Hogg, A.J. and Pritchard, D., “The effects of hydraulic resistance on
dam-break and other shallow inertial flows”, J. Fluid Mech., 501(2004), pp. 179-202.
(34) Liu, T.P., “Hyperbolic conservation laws with relaxation”, Commun.
Math. Phys., 108(1) (1987), pp. 153-175.
(35) Berthon, C., Chalons, C. and Turpault, R., “Asymptotic-preserving
Godunov-type numerical schemes for hyperbolic systems with stiffand non-stiff relaxation terms”, Numer. Methods Partial Differ.
Equat., 29(4) (2013), pp. 1149-1172.
(36) Yoshioka, H., Unami, K. and Fujihara, M., “Comparative numericalanalysis on momentum flux evaluation schemes for shallow waterflows in open channel networks”, J. Rainwater Catch. Syst. (inpress)(37) Ashida, K. and Kawai, S., “Hydraulic behaviors of dividing flow
around junctions”, Disaster Prevention Research Institute Annuals,
55 (1979), pp. 491-506. (in Japanese)(38) Koseki, H., Hoshino, T. and Yasuda, H., “Study on distribution ratioin bifurcation of response to down stream boundary condition”, J.
JSCE, Ser. B1, 69(4) (2013), pp. I_601-I_606. (in Japanese)(39) Goutal, N. and Maurel, F., Proc. the 2nd Workshop on Dam-break
Wave Simulation, EDF-DER Report, HE-43/97/016/B, 1997.
(40) Vázquez-Cendón, M.L., “Improved treatment of source terms in
upwind schemes for the shallow water equations in channels withirregular geometry”, J. Comput. Phys., 148(2) (1999), pp. 11-34.
(41) Zoppou, C. and Roberts, S. “Explicit schemes for dam-break flows”,
J. Hydraul. Eng., 129(1) (2003), pp.11-34.

10

Copyright © 2013 by JSFM