ACKNOWLEDGEMENTS

Foremost, I bow my head with deep gratitude to Almighty Allah Who bestowed upon me
the courage and steadfastness to accomplish my under graduation. He is the most powerful,
compassionate, kind and merciful. I offer my most humble words of thanks to the Holy
Prophet Muhammad (peace be upon him) who is forever a torch of guidance for humanity.
I would like to express my deep and sincere gratitude to my supervisor Dr. habil. Shamsul
Qamar Associate Professor at Department of Mathematics, COMSATS Institute of Infor-
mation Technology, Islamabad, Pakistan. He is also working as visiting professor at Max
Planck Institute for Dynamics of Complex Technical Systems Magdeburg, Germany. I am
grateful to him for his devotion, keen interest in understanding, encouragement, enthusi-
asm, personal guidance, immense knowledge and motivation. His guidance helped me in all
the time of research and writing of this thesis. I am deeply grateful to him for his detailed
and constructive comments. His logical way of thinking has been of great value for me
throughout this work. I could not have imagined having a better advisor and mentor for
my thesis. Indeed, it is a matter of pride and privilege for me to be his graduate student.
Along with my worthy supervisor, I wish to acknowledge Mr. Saqib Zia for being a great
helping hand along with his invaluable and intellectual guidance, beneficial remarks, con-
structive criticism, encouragement, insightful comments, and hard questions.
Most importantly I am grateful to COMSATS Institute of Information Technology, Islam-
abad and the worthy Rector Dr. S. M. Junaid Zaidi for providing excellent opportunities.
I am thankful to all my fellow mates for extending a hand full of help whenever required
also for the stimulating discussions and for the sleepless nights we were working together
before deadlines.
Last but not the least, I would like to thank my family, especially my parents for being
my loving spiritual and motivational support throughout my life. Without their encour-
agement and understanding it would have been impossible for me to finish this work.

Nadeem Ahmed
CIIT/FA11-RMT-006/ISB

viii
 ABSTRACT

Application of Space-Time CESE Method for

Solving Two-Phase Shallow Flow Model

with Variable Bottom Topography

This thesis is related to numerical simulation of (SPSF)and (TPSF) models. A space time
CESE scheme is suggested to solve these models. In the absence of bottom topography the
SPSF form a homogenous hyperbolic system while TPSF model governs a non-homogenous
system of conservation laws. Therefore, the proposed scheme is derived for both types of
systems. Numerous applications of CESE method in different areas reveal the method’s
generality, feasibility and effectiveness. A number of numerical test problems are presented
in this thesis. The results of proposed scheme are compared with those obtained from
the central scheme. The comparison shows better performance of the CESE scheme in
simulating the single and two-phase flows.

ix
Contents

1 Introduction 1
1.1 Applications of the Flow Model . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.3 New Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.1 Law of Conservation of Mass . . . . . . . . . . . . . . . . . . . . . . 6
1.4.2 Law of Conservation of Momentum . . . . . . . . . . . . . . . . . . 6
1.4.3 Shallow Water Equations . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.4 Multiphase Flow . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.4.5 The Scalar Conservation Law . . . . . . . . . . . . . . . . . . . . . 7
1.4.6 Hyperbolic Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4.7 The Riemann Problem . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4.8 Weak Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.5 Thesis Layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2 Derivation of Single and Two-Phase Shallow Flow models 13

2.1 Derivation of the Equation of Continuity . . . . . . . . . . . . . . . . . . . 13
2.2 Derivation of Conservation of Momentum . . . . . . . . . . . . . . . . . . . 14
2.3 Derivation of Single-Phase Shallow Flow Model . . . . . . . . . . . . . . . 19
2.3.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Derivation of Two-Phase Shallow Flow Model . . . . . . . . . . . . . . . . 27
2.4.1 Steady States . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

x
 2.5 The Two-Dimensional Two-Phase Shallow Granular Flow Model . . . . . . 34

3 The CESE Scheme 36

3.1 Numerical Scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.2 The CESE Method for Single-Phase Shallow Flows . . . . . . . . . . . . . 36
3.3 The CESE Method for Two-Phase Shallow Flows . . . . . . . . . . . . . . 41

4 Case Studies 46
4.1 Single-Phase Shallow Flows . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.1.1 Test Problem 1: Rarefaction into Vacuum . . . . . . . . . . . . . 46
4.1.2 Test Problem 2: Formation of Dry bed . . . . . . . . . . . . . . . 46
4.1.3 Test Problem 3: Test with Variable Bottom Topography . . . . . . 47
4.2 Test Problems of Two-Dimensional Single-Phase Shallow Flows . . . . . . 47
4.2.1 Test Problem 4: Two Separated Fluids . . . . . . . . . . . . . . . 47
4.2.2 Test Problem 5: Rotor Like Problem . . . . . . . . . . . . . . . . 51
4.3 Two-Phase Shallow Flows . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.3.1 Test Problem 6: Simple Riemann Problem . . . . . . . . . . . . . 54
4.3.2 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.3.3 Test Problem 7: Rarefaction of the Fluid Constituent into Vacuum 54
4.3.4 Test Problem 8: Spreading of a Granular Mass . . . . . . . . . . . 56
4.3.5 Test Problem 9: Generation of Dry Bed . . . . . . . . . . . . . . 57
4.3.6 Test Problem 10: Interface Propagation . . . . . . . . . . . . . . . 58
4.4 Test Problem with Variable Bottom Topography . . . . . . . . . . . . . . . 58
4.4.1 Perturbation of a Steady State at Rest . . . . . . . . . . . . . . . . 58
4.4.2 Two-Dimensional Two-Phase Problem . . . . . . . . . . . . . . . . 59

5 Conclusion 70

Bibiliography 75

xi
List of Figures

1.1 Classification of two-phase flows. . . . . . . . . . . . . . . . . . . . . . . . 8

1.2 Flow across two points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.1 Space-time staggered grid near SE (j, n) . . . . . . . . . . . . . . . . . . . 39

4.1 Rarefaction into vacuum for SPSF. . . . . . . . . . . . . . . . . . . . . . . 48

4.2 Dry bed formation for SPSF at t=1. . . . . . . . . . . . . . . . . . . . . . 49
4.3 Problem 3: Results of variable topography at t=8 and t=12. . . . . . . . . 50
4.4 Problem 4: Two separated fluids results on 100 × 100 grid points at t=0.15. 52
4.5 Problem 5: Rotor like problem results on 100 × 100 grid points t = 0.15. . 53
4.6 : Error Analysis for height h and solid volume fraction φ. . . . . . . . . . . 55
4.7 Test Problem 5: Comparison of results at t = 0.5 over 200 cells. . . . . . . 60
4.8 Test Problem 7: Comparison of results over 100 cells. . . . . . . . . . . . . 61
4.9 Test Problem 8: Comparison of results for granular mass spreading. . . . . 62
4.10 Test Problem 9 (i): Comparison of results for generation of dry bed problem. 63
4.11 Test Problem 9 (ii), Comparison of results over domain [-5,5] at time t=1 . 64
4.12 Test Problem 9 (iii), Results of for height, momenta and velocities. . . . . 65
4.13 Problem 10 : Results for Interface evolution, total height and velocities. . . 66
4.14 Problem 11 : Perturbation of steady-state at rest (h̃ = φ̃ = 10−3 ). . . . . . 68
4.15 Problem 10: Results on 100 × 100 grid points at time t = 0.5. . . . . . . . 69

xii
List of Tables

4.1 L1 -errors for phi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.2 L1 -error velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

xiii
LIST OF ABBREVIATIONS

CESE Conservation Element and Solution Element

CEs Conservation Elements
CFD Computational Fluid Dynamics
CFL Courant-Friedrich-Lewy
FDM Finite Difference Method
FEM Finite Element Method
FVM Finite Volume Method
SEs Solution Elements
SPSF Single-Phase Shallow Flow
SWEs Shallow Water Equations
SWWEs Shallow Water Wave Equations
TPSF Two-Phase Shallow Flow

xiv
LIST OF NOTATIONS

b Bottom Topography
ds Element Surface Area
dσ Area
F Flux
g Gravitational Acceleration
h Flow Height
e
H Characteristic Depth
e
L Horizontal Scale
n unit outward normal vector
S Surface
v Fluid Velocity
V Volume
W Conserved Quantity
η Interface Between Each Layer
ξ Intermediate Surface Coordinate
φ Solid Volume Fraction
ρ Density
τ Stress Tensor

xv
Chapter 1

Introduction

Many types of flow, not necessarily involving water as a fluid, can be characterized as
shallow-water flows. The general characteristic of such flows is that the vertical dimension
is much smaller than any typical horizontal scale and this is true in many everyday situ-
ations. These equations arise from the basic equations of fluid mechanics for an inviscid
and incompressible fluid. The basic equations governing fluid motion are based on the
principles of conservation of mass, momentum and energy.
Gravitational geophysical flows typically involve both solid granular material and intersti-
tial fluid. Our main task is to numerically approximate the two-phase shallow flow (TPSF)
model. The model describes flow of solid and fluid mixtures, consisting of momentum and
mass equations for both the phases. It corresponds to the well-known Pitman-Le two-fluid
system coupled together by momentum exchange terms.
The relative simplicity of the model is the basic reason of attraction for many mathemati-
cians to study such flows. This is not a historic review, for more historical information see
e.g book by Dronkers [1]. Some names should however be mentioned because of terminol-
ogy. A lot of interesting work on tidal flows was done by Laplace, known as the Laplacian
tidal equations which is a specialized form of the shallow water equations [2]. In the French
scientific community, the SWEs are commonly referred to as the Saint-Venant equations,
although it appears that Saint-Venant derived only the 1-D version [1]. Here we will use
the term shallow water equations (SWEs) throughout the thesis.

1
The numerical solution of the SWEs was one of the early applications of digital computers
when these became available in the late 1940’s. Simulations were done by Charney et al.
[3] for atmospheric flows. A considerable development has been taken place since its intro-
duction and it can be stated that at least 2-D shallow water equation are well understood.
Shallow water flows are characterized by flow regions with a free surface, an impermeable
bottom topography such as a sea floor and the horizontal velocity that dominates the flow
field. Such flows may be generated by gravitational forces associated with sloped beds,
wind that creates stress on the free surface of the fluid or an applied pressure gradient such
as a tidal influence, a disturbance in the impermeable bottom say from an earthquake or
a combination of two ar all three.

1.1 Applications of the Flow Model

The two-phase SWEs are widely used in oceans and in hydraulic engineering to model flow
in rivers, reservoirs and coastal areas. Similarly, shallow water wave equations have been
used for many applications and some of them are listed below.

• Storm Surges
Storm surges happens due to large wind velocities, pressure and high tides associated
with large storms, especially tropical cyclones. It is important to be able to model
storm surges to identify areas of high risk and areas that may be protected by artificial
barriers such as break walls. Examples of modeling storm surges can be found in [4, 5]

• Solute Transport
Prediction of solute transport concentrations in various fluid flows is important for
many environmental and industrial applications. When developing ecological models
of wetlands and estuaries information regarding the movement of pollutants and other
form of macro solutes which are often food sources is invaluable. See [6] for details
on the use of the shallow water wave equations (SWWEs) for modeling the transport

2
 of the passive pollutants. Other examples can be found in [7, 8]

• Tsunami Propagation
Tsunami poses a large risk to coastal communities all over the world. This became
agonizingly evident in the aftermath of the so-called boxing-day tsunami of 2004.
For these reasons there is an increased focus on tsunami hazard mitigation (tsunami
detection, forecasting and emergency preparedness) to limit loss of life and economic
fallout resulting from such an event [9]. The SWEs are commonly used to model
tsunami propagation [10, 11].

• Dam Breaks
Water impoundment resulting from dam breaks poses a potential risk to life and
infrastructure. Such events are often characterized by rapidly varying flows and
shock waves. In [12] the authors have used the shallow flow model to simulate the
sudden collapse of supply reservoirs in urban areas.

• River Flows and Inundation

Flow in rivers over flood plains is very accurately represented by shallow water flow
model. In [13], such flows were simulated to determine flow patterns and river heights
during flooding. In [14], the shallow water flow theory was used to model steady flow
in rivers.

• Ecological Models
Macro particles found in rivers and oceans are often food sources for small marine
animals. The concentration of these particles can thereby have an effect on the pop-
ulation size and distributions of these organisms and in turn on the ecosystem as a
whole. Consequently, shallow water modeling can be used to simulate the movement
of these particles and organisms and used to determine and predict important eco-
logical markers. In [15] the authors have used shallow water model coupled with a
dispersion model to predict the distribution of prawn larvae in the Spencer Gulf.

3
1.2 Literature Review
Savage and Hutter [16] initiated the work on shallow water equations in one dimension.
Various authors extended this approach to two-space dimensions and flows having complex
bottom topography.

The proposed two-phase shallow flows (TPSF) model was first discussed by Pelanti and
co-workers in [17], where they proved that the model is hyperbolic for sufficiently small
difference of velocity between the phases (fluid and solid). By this assumption of hyper-
bolicity, the current TPSF model has great applications to geophysical flows,e.g. debris
flows and land slides where the drag forces between the phases rapidly move fluid and
solid constitutes towards the state of kinematic equilibrium [18]. Initially, the TPSF model
was numerically solved by finite volume schemes. This scheme is reliant upon Roe-type
Riemann solver. The scheme may generate unphysical negative values of height of flow as
well as of phase volume fraction and this is one of the disadvantage of the scheme. Positiv-
ity preserving numerical schemes can handle efficiently wet/dry transition regions in free
surface flow models.

There are many Riemann solvers which preserve positivity for single-phase shallow flow
(SPSF) model. Some of them are the exact Riemann solvers, for instance well known
Riemann solvers (HLL, HLLC) [19, 20]. Other Relaxation based solver include the recent
scheme given in [21]. This scheme [21] have the efficiency to deal with vacuum states. One
of the methods proposed by Castro and co-workers was modified Roe method (MRoe) This
method does not preserve positivity rigorously to avoid negative water depths, it requires
restriction on the CFL number.

Unluckily, the above discussed positivity preserving solvers can not be extended to the
TPSF model presented in this thesis. The nonlinearly coupled hyperbolic system of the
model has non-conservative flux terms. Also the solution of the model for Riemann struc-
ture is quite complex. Major problems for the current TPSF model include the implicit
expression for eigenvalues, the non-availability of Riemann invariants and the deficiency of

4
knowledge about the exact Riemann solution structure.

A relaxation strategy based Riemann solver was developed by Berthon and co-author
[21] for shallow flow models. In [21] they preserved the positivity of their method for
non-conservative variables. The approach of positivity in [21] is inspired by the Riemann
Invariants developed for SPSF model. However, this scheme can not be extended for the
proposed TPSF model. We do not know about the Riemann invariants of the model and
it is no more analogous to the model given in [21]. Therefore, positivity of the current
TPSF model is not ensured by the argument discussed in [21]. In [22], the authors have
rigorously preserved, the positivity of their scheme for SPSF model but the scheme does
not ensure positivity for the current TPSF model in general and an appropriate condition
is needed on the CFL number.

1.3 New Contribution

In this thesis, the space-time CESE method of Chang [23] is extended to solve SPSF and
TPSF models. It is based on a different concept, differs from the familiar methods e.g.
the FDM, FVM and FEM. It has many unique features, e.g space and time are treated
in a unified fashion, no Riemann solver is needed in calculating interfacial fluxes and
systems of conservation laws are treated in exactly the same way as single conservation laws.
Numerous applications of CESE method in different areas reveal the method’s generality,
feasibility and effectiveness. These areas include unsteady flows problems [23, 24, 25],
vortex dynamics in aeroacoustics [26], problems of diffusion [27], flows (viscous) [28], jets
(supersonic) [29], axisymmetric and inviscid flows [30] and magnetohydrodynamics [31].
We follow a similar procedure presented in [16, 36] to make the current scheme well balanced
in the case of variable bottom topography. In the end, CESE method was compared with
staggered central scheme [33] and results were compared.

5
1.4 Preliminaries
This section introduces some basic concepts and definition to understand the entire work
presented in this thesis.

1.4.1 Law of Conservation of Mass

The law states that mass can neither be created nor destroyed although the entities related
to it may change their form. The law states that for any system, the mass of the system
must remain constant over time, as system mass will not change if there is no addition or
removal.

1.4.2 Law of Conservation of Momentum

This law states that the total momentum of an isolated system before collision is always
equal to the total momentum after collision means when some bodies constituting an
isolated system act upon one another, the total momentum of the system remains constant

1.4.3 Shallow Water Equations

The shallow water equations are a set of hyperbolic equations for the full free-surface
problem. In the uni-dimensional form these equations are also known as Saint Venant
equations. These equations are obtained from depth integrating Navier-Stokes equations
where the horizontal length scale is much lager than the vertical length scale.
e is the characteristic depth and L
If H e horizontal scale of motion, then for shallow water
equations
e
H
δ= ≪ 1. (1.1)
e
L

1.4.4 Multiphase Flow

Multi-phase flow is a fundamental mechanism of CFD. It is a process in which two or more
fluids flow are involved in a region, specifically two in this case. We get encountered to

6
the two-phase flow in various scientific and technical fields, e.g. environmental research,
chemical and processing engineering installations and nuclear operations. In daily life,
boilers and geysers are the practical examples of two-phase flows. According to the phase
materials, two phase flows can be subdivided into three categories: gas-liquid flows (bub-
bly flows, separated flows, gas-droplet flows), solid-gas flows (gas-particle flows, fluidized
beds), liquid-solid flows (slurry flows, sediment transport) as illustrated in Fig. 1.1.

1.4.5 The Scalar Conservation Law

Partial differential equation of the form

Wt + F (W )x = 0, (1.2)

is called a scalar conservation law in one space dimension. Here, W = W (x, t) is known as
the conserved quantity, while F is the flux. The variable t denotes time, while x is the one
dimensional space variable.

Equation (1.2) often describes transport phenomena. On integrating the equation (1.2)
over a given interval [a, b], we obtain

Zb Zb Zb
d
W (x, t)dx = Wt (x, t)dx = − F (W (x, t))x dx, (1.3)
dt
a a a

= F (W (a, t)) − F (W (b, t)), (1.4)

= [Inflow at a] - [outflow at b]. (1.5)

In other words, if the inflow is equal to the outflow then the quantity W is neither created
nor destroyed. The total amount of W contained inside any given interval [a,b] can change

7
Figure 1.1: Classification of two-phase flows.

8
 u

a b ξ x

Figure 1.2: Flow across two points.

only due to the flow of W across the boundary points as shown in Fig. 1.2.

1.4.6 Hyperbolic Systems

The models studied in this thesis are the hyperbolic systems of partial differential equations
(PDEs). Consider the one-dimensional conservation law as

∂W ∂F(W)
+ = 0, (1.6)
∂t ∂x

where,
   
W F1
 1  
   
 W2   F2 
W= 
 ..  , F(W) = 
 ..
.
 (1.7)
 .   . 
   
Wm Fm

Here, F(W) and W are vectors of fluxes and conserved variables respectively. Moreover
m represents the number of equations in the system. The chain rule allows us to reduce
the system (1.6) into quasi-linear form

∂W ∂F(W)
+ F′ (W) = 0. (1.8)
∂t ∂x

9
The Jacobian matrix of flux functions are given as
 
∂F1 ∂F1
,··· , ∂W
∂F(W)   ∂W1 m
.. .. .. 
A(W) = = . . .  (1.9)
∂W  
∂Fm ∂Fm
∂W1
,··· , ∂Wm

The entries of matrix A(W) are the partial derivatives of the components of the flux vector
F with respect to the components of the vector of conserved variables W.
The eigenvalues λi , i = 1, 2, · · · , m, of matrix A are the solutions of characteristic polyno-
mials

|A − λI| = det(A − I) = 0 , j = 1, 2, · · · , m. (1.10)

Here, we used usual notation of identity matrix as I.

The system (1.6) is said to be hyperbolic at a point (x, t) if A has m real eigenvalues
λ1 , · · · , λm with a corresponding set of m linearly independent right eigenvectors. The
system is said to be strictly hyperbolic if the eigenvalues λi are all distinct.

1.4.7 The Riemann Problem

The solution to Riemann problems of shallow water equations are investigated in this
work. In this problem the system has discontinuous initial data given by two constant
states i.e. a left W l and a right W r , separated by a plane across which at least one of the
fields experiences an instantaneous jump, i.e. for at least one i, Wil 6= Wir . The resulting
problem is known as the Riemann problem. The study of a general Riemann problem
allows one to investigate the behavior of weak solutions.

1.4.8 Weak Solution

Let W : R 7→ R be a smooth function. A measurable function W = W (t, x), defined on an
open set Ω ⊆ R+ × R and with values in R, is a weak solution of the scalar conservation

10
laws
Wt + F (Wx ) = 0, (1.11)

if, for every C 1 function φ : Ω 7→ R with compact support, one has

Z Z
{W φt + F (W )φx }dxdt = 0. (1.12)
Ω

Here, no continuity assumption is made on W . One only need W and F (W ) to be locally

integrable in Ω. Another form of weak formulation is called Oleinik’s integral formula
which is given as,

I
W dx − F (W )dt = 0. (1.13)
∂Ω

Here ∂Ω is the boundary of the domain.

1.5 Thesis Layout

The remaining thesis is arranged in the following manner.

Chapter 2
In this chapter, we have derived conservation of mass and conservation of momentum. This
chapter, introduces one and two dimensional SPSF models. The models are extended to
one and two dimensional two phase shallow flows (TPSF) and their derivation is included.

Chapter 3
In this chapter, we give a brief introduction to the CESE for SPSF and TPSF models and
derive the suggested CESE scheme for solving these models.

Chapter 4
In chapter 4, several test problems are carried out for all the two types models. At the

11
end, comparison of the suggested CESE scheme and central scheme is carried out.

Chapter 5
This chapter presents thesis conclusion and further recommendations.

12
Chapter 2

Derivation of Single and Two-Phase

Shallow Flow models

In this chapter we will derive SPSF and TPSF models. First we will derive equation of
continuity and momentum. Afterwards, we will go the derivation of models.

2.1 Derivation of the Equation of Continuity

Consider a fluid with local velocity v and local density ρ. Suppose a control volume V has
a boundary S then the total mass is given by,

Z
M= ρdV (2.1)

Rate of change of this mass is given by,

Z
∂M ∂ρ
= dV (2.2)
∂t ∂t

This happens due to stuff flowing across the boundary therefore,

13
 Z
∂M
= ρv.dS (2.3)
∂t

Now by applying Green’s theorem, we get

Z
∂M
=− ∇.(ρv)dV (2.4)
∂t

Now Eq. (2.2) and (2.4) are compared to get integrands to be equal. From here we get the
local conservation of mass.

∂ρ
+ ∇.(ρv) = 0 (2.5)
∂t

The above equation is also called a continuity equation.

2.2 Derivation of Conservation of Momentum

The equation representing the law of conservation of momentum is deduced from the New-
ton’s second law of motion:

X dv
F=m = ma, (2.6)
dt

since,

m = ρdV,

m = ρdxdydz, as dV = dxdydz.

14
Now,
d dv
(mv) = m ,
dt dt
dv
= ρdxdydz .
dt

The right hand side of above equation can be written as,

dv d
ρdxdydz = ρdxdydz (uî + v ĵ + w k̂),
dt dt
d
= ρdxdydz (u v w),
dt
X
= (Fx Fy Fz ).

where Fx , Fy , Fz are components of vector F. In component form, the above equation can
be written as,
du X
ρdxdydz = Fx , (2.7)
dt
dv X
ρdxdydz = Fy , (2.8)
dt
dw X
ρdxdydz = Fz . (2.9)
dt

Now,

v = v(x, y, z, t),
∂v ∂v ∂v ∂v
dv = dx + dy + dz + dt,
∂x ∂y ∂z ∂t
dv ∂v dx ∂v dy ∂v dz ∂v
= + + + ,
dt ∂x dt ∂y dt ∂x dt ∂t
∂v ∂v ∂v ∂v
=u +v +w + ,
∂x ∂y ∂z ∂t
∂ ∂ ∂ ∂v
= (u v w).( )v + .
∂x ∂y ∂z ∂t

15
In compact form this can be written as,

dv ∂v
= (v.∇)v + . (2.10)
dt ∂t

From above equation, change in momentum in x-direction is given by,

du ∂u
= (v.∇)u + , (2.11)
dt ∂t

dv ∂v
= (v.∇)v + , (2.12)
dt ∂t

dw ∂w
= (v.∇)w + . (2.13)
dt ∂t

Forces which are acting are body forces, surface forces and pressure i.e.

X X
Fx = (Bx + Px + Sx ), (2.14)
X X
Fy = (By + Py + Sy ), (2.15)
X X
Fz = (Bz + Pz + Sz ). (2.16)

If bx is component of body force per unit volume,

Bx = ρdxdydzbx (2.17)

By = ρdxdydzby , (2.18)

Bz = ρdxdydzbz . (2.19)

In compact form, this can be written as,

B = ρdxdydz(bx by bz ). (2.20)

16
On combining three components we have,

d dv
(uî + v ĵ + w k̂) = = ρf − ∇P + ∇.S, (2.21)
dt dt

where
 
S S S
 xx yx zx 
 
S =  Sxy Syy Szy .
 
Sxz Syz Szz

Now,

∂v
ρ[ + (v.∇)v] = ρf − ∇P + ∇S. (2.22)
∂t

In component form, above equation can be written as,

∂u ∂u ∂u ∂u ∂P ∂Sxx ∂Syx ∂Szx

ρ[ +u +v + w ] = ρbx − + + + , (2.23)
∂t ∂x ∂y ∂z ∂x ∂x ∂y ∂z
∂v ∂v ∂v ∂v ∂P ∂Sxy ∂Syy ∂Szy
ρ[ +u +v + w ] = ρby − + + + , (2.24)
∂t ∂x ∂y ∂z ∂y ∂x ∂y ∂z
∂w ∂w ∂w ∂w ∂P ∂Szx ∂Syz ∂Szz
ρ[ +u +v +w ] = ρbz − + + + . (2.25)
∂t ∂x ∂y ∂z ∂z ∂x ∂y ∂z

Now,

∂2S ∂u
= 2µ ,
∂x2 ∂x
∂2S ∂u ∂v ∂2S
= µ( + )= ,
∂x∂y ∂y ∂x ∂y∂x
∂2S ∂u ∂w ∂2S
= µ( + )= ,
∂x∂z ∂z ∂x ∂z∂x
∂2S ∂v ∂w ∂2S
= µ( + )= .
∂y∂z ∂z ∂y ∂z∂y

17
where µ is viscosity. From momentum equation, x-component is given by,

∂u ∂u ∂u ∂u ∂P ∂ ∂u ∂ ∂u ∂v ∂ ∂w ∂u
ρ[ +u +v + w ] = ρbx − + (2µ ) + [µ( + )] + [µ( + )].
∂t ∂x ∂y ∂z ∂x ∂x ∂x ∂y ∂y ∂x ∂z ∂x ∂z

Similarly y-component is,

∂v ∂v ∂v ∂v ∂P ∂ ∂u ∂v ∂ ∂v ∂ ∂v ∂w
ρ[ +u +v + w ] = ρby − + [µ( + )] + (2µ ) + [µ( + )],
∂t ∂x ∂y ∂z ∂y ∂x ∂y ∂x ∂y ∂y ∂z ∂z ∂y

and z-component becomes,

∂w ∂w ∂w ∂w ∂P ∂ ∂w ∂u ∂ ∂v ∂w ∂ ∂w
ρ[ +u +v +w ] = ρbz − + [µ( + )] + [µ( + )] + (2µ ).
∂t ∂x ∂y ∂z ∂z ∂x ∂x ∂z ∂y ∂z ∂y ∂z ∂z

From x-component we have,

∂u ∂u ∂u ∂u ∂P ∂2u ∂2u ∂2v ∂2w ∂2u

ρ[ +u +v + w ] = ρbx − + 2µ 2 + µ 2 + µ +µ +µ 2,
∂t ∂x ∂y ∂z ∂x ∂x ∂y ∂y∂x ∂z∂x ∂z
∂P ∂2u ∂2u ∂2u ∂2u ∂2v ∂2w
= ρbx − + µ( 2 + 2 + 2 ) + µ( 2 + + ),
∂x ∂x ∂y ∂z ∂x ∂y∂x ∂z∂x
∂P ∂ ∂u ∂v ∂w
= ρbx − + µ∇2 u + µ ( + + )
∂x ∂x ∂x ∂y ∂z
∂P ∂
= ρbx − + µ∇2 u + µ (∇.v).
∂x ∂x

For incompressible viscous fluid,

∇.v = 0. (2.26)

Thus we have,
∂u ∂u ∂u ∂u ∂P
ρ[ +u +v + w ] = ρbx − + µ∇2 u, (2.27)
∂t ∂x ∂y ∂z ∂x

18
which is x-component of momentum equation for incompressible viscous (Newton) fluids.
Similarly y-component is,

∂v ∂v ∂v ∂v ∂P
ρ[ +u +v + w ] = ρby − + µ∇2 v, (2.28)
∂t ∂x ∂y ∂z ∂y

and z-component is,

∂w ∂w ∂w ∂w ∂P
ρ[ +u +v +w ] = ρbz − + µ∇2 w. (2.29)
∂t ∂x ∂y ∂z ∂z

2.3 Derivation of Single-Phase Shallow Flow Model

We start our discussion with the construction of a depth averaged model using the three
dimensional, inviscid Navier-Stokes equations [34]:


 ∂ρ
 + ∇.uρ = 0
∂t (2.30)
 ∂ (ρu) + ∇.(ρuu) = −∇P + ∇.τ − ρgẑ

∂t

where u ≡ (u, v, w)T is the velocity in the three direction, P is the pressure, g is the accel-
eration due to gravity and τ is the stress tensor. The boundary conditions are,

19
 ∂η ∂η ∂η
w= +u +υ and P = PA (x, y, t) at z = η,
∂t ∂x ∂y
∂2τ ∂ 2 τ ∂η ∂ 2 τ ∂η ∂2τ
=− 2 − + at z = η,
∂s∂x ∂x ∂x ∂x∂y ∂y ∂x∂z
∂2τ ∂ 2 τ ∂η ∂ 2 τ ∂η ∂2τ
=− − + at z = η,
∂s∂y ∂y∂x ∂x ∂y 2 ∂y ∂y∂z
∂2τ ∂ 2 τ ∂b ∂ 2 τ ∂b ∂2τ
= 2 + − at z = η,
∂b∂x ∂x ∂x ∂x∂y ∂y ∂x∂z
∂2τ ∂ 2 τ ∂b ∂ 2 τ ∂b ∂2τ
= + − at z = η,
∂b∂y ∂x∂y ∂x ∂y 2 ∂y ∂y∂z
∂b ∂b
w =u +υ at z = b,
∂x ∂y

representing a wide range of permissible geophysical flows with a top free surface η subject
to stresses τsx , τsy and an impermeable bottom boundary whose elevation is determined
from a given function b(x, y, t) and subject to stresses τbx , τby . In addition to the stresses
at the top surface, we also will assume that an atmospheric pressure at the top surface
PA (x, y, t) is also given. If we also make the simplifying assumption that ρ is constant, the
advection of density in (2.1) simplifies to

∇.u = 0. (2.31)

This allows us rewrite the momentum conservation equations as

∂u 1
+ ∇.(uu) = [−∇P + ∇.τ ] − gẑ. (2.32)
∂t ρ

Next, we state the Liebniz Theorem [35] which will be used in the derivation of the model.

Liebniz Theorem

Z b(y,t) Z b(y,t)
∂ ∂f ∂a ∂b
f (x, y, t)dx = dx − f (a, y, t) + f (b, y, t) . (2.33)
∂t a(y,t) a(y,t) ∂t ∂t ∂t

20
2.3.1 Derivation

Applying Leibniz theorem and integrating the continuity equation through the depth leads
to

Z η Z η Z η
∂u ∂v ∂w ∂ ∂ ∂η ∂η ∂b ∂b
+ + dz = udz + υdz + (w − u − v )|z=η − (w − u − υ )|z=b
b ∂x ∂y ∂z ∂x b ∂y b ∂x ∂y ∂x ∂y
∂ ∂ ∂η
= (hu) + (hυ) + =⇒
∂x ∂y ∂t
∂h ∂ ∂
= + (hu) + (hυ),
∂t ∂x ∂y

where we have used η = h + b and defined

Z η
1
f≡ f dz,
h b

representing the average of the quantity f through the depth. For the momentum equa-
tions, the advective terms in the x-direction are

Z η Z Z η
η Z η
∂u ∂ 2 ∂ ∂ ∂ ∂ 2 ∂
[ + (u ) + (uυ) + (uw)]dz = udz + u dz + uvdz
b ∂t ∂x ∂y ∂z ∂t b ∂x b ∂y b
∂η ∂η ∂η
+ [u(w − u −v − )]z=η
∂x ∂y ∂t
∂b ∂b
− [u(w − u − υ )]z=b
∂x ∂y
∂ ∂ ∂
= (hu) + (hu2 ) + (huv).
∂t ∂x ∂y

For the y-direction

21
 Z η Z η Z η Z η
∂v ∂ ∂ 2 ∂ ∂ ∂ ∂
[ + (uυ) + (υ ) + (υw)]dz = υdz + uυdz + υ 2 dz
b ∂t ∂x ∂y ∂z ∂t b ∂x b ∂y b
∂η ∂η ∂η
+ [v(w − u −v − )]z=η
∂x ∂y ∂t
∂b ∂b
− [υ(w − u − υ )]z=b
∂x ∂y
∂ ∂ ∂
= (hυ) + (huυ) + (hυ 2 ).
∂t ∂x ∂y

Finally for the z-momentum

Z η Z η Z η Z η
∂w ∂ ∂ ∂ ∂ ∂ ∂
[ + (uw) + (υw) + (w 2 )]dz = wdz + uwdz + υwdz
b ∂t ∂x ∂y ∂z ∂t b ∂x b ∂y b
∂η ∂η ∂η
+ [w(w − u −v − )]z=η
∂x ∂y ∂t
∂b ∂b
− [w(w − u − υ )]z=b
∂x ∂y
∂ ∂ ∂
= (hw) + (huw) + (hυw).
∂t ∂x ∂y

Turning to the right side of the momentum conservation equation (2.31),

1
[−∇P + ∇.τ ].
ρ

Let us first consider the pressure term P . Assume that P is of the form

P (x, y, z, t) = PA (x, y, t) + ρg(η − z) + p(x, y, z, t),

where again PA is the pressure at the surface of the flow, the second term is due to the
hydrostatic pressure and p is the deviation from hydrostatic pressure. Inserting these into

22
the gradient leads to

 
∂ ∂η
(PA ) + ρg ∂x
 ∂x 
 ∂ ∂η 
∇P =  (PA ) + ρg ∂y  + ∇p. (2.34)
 ∂y 
−ρg

Integrating ∇P through the depth then gives

 
∂ ∂η
Z (PA ) + ρg ∂x
η  ∂x  Z η
 
∇P dz = h  ∂
(PA ) + ρg ∂η + ∇pdz. (2.35)
b  ∂y ∂y  b
−ρg

The stress tensor τ can also be integrated in each direction where in the x-direction we
have

Z η Z η 2 Z η 2
∂τxx ∂τxy ∂τxz ∂ ∂ τ ∂ ∂ τ ∂η ∂ 2 τ ∂η ∂ 2 τ ∂2τ
( + + )dz = dz + dz − ( + − )
b ∂x ∂y ∂z ∂x b ∂x2 ∂y b ∂x∂y ∂x ∂x2 ∂y ∂x∂y ∂x∂z
∂b ∂ 2 τ ∂b ∂ 2 τ ∂2τ
+( + − ),
∂x ∂x2 ∂y ∂x∂y ∂x∂z
Z η 2 Z η 2
∂ ∂ τ ∂ ∂ τ ∂2τ ∂2τ
= dz + dz + + ,
∂x b ∂x2 ∂y b ∂x∂y ∂s∂x ∂b∂x

and in the y-direction we have

Z η Z η 2 Z η 2
∂τyx ∂τyy ∂τyz ∂ ∂ τ ∂ ∂ τ ∂η ∂ 2 τ ∂η ∂ 2 τ ∂2τ
( + + )dz = dz + dz − ( + − )
b ∂x ∂y ∂z ∂x b ∂y∂x ∂y b ∂y 2 ∂x ∂y∂x ∂y ∂y 2 ∂y∂z
∂b ∂ 2 τ ∂b ∂ 2 τ ∂2τ
+( + − ),
∂x ∂y∂x ∂y ∂y 2 ∂y∂z
Z η 2 Z η 2
∂ ∂ τ ∂ ∂ τ ∂2τ ∂2τ
= dz + dz + + .
∂x b ∂y∂x ∂y b ∂y 2 ∂s∂y ∂b∂y

23
Finally, in the z-direction we have
Z η Z η 2 Z η 2
∂τzx ∂τzy ∂τzz ∂ ∂ τ ∂ ∂ τ ∂η ∂ 2 τ ∂η ∂ 2 τ ∂2τ
( + + )dz = dz + dz − ( + − 2 )z=η
b ∂x ∂y ∂z ∂x b ∂z∂x ∂y b ∂z∂y ∂x ∂z∂x ∂y ∂z∂y ∂z
2 2 2
∂b ∂ τ ∂b ∂ τ ∂ τ
+( + − 2 )z=b ,
∂x ∂z∂x ∂y ∂z∂y ∂z
∂ ∂ ∂η ∂ 2 τ ∂η ∂ 2 τ ∂2τ
= (hτzx ) + (hτzy ) − ( + − 2 )z=η
∂x ∂y ∂x ∂z∂x ∂y ∂z∂y ∂z
2 2 2
∂b ∂ τ ∂b ∂ τ ∂ τ
+ + )z=b ,
∂x ∂z∂x ∂y ∂z∂y ∂z 2

where τsx ,τsy ,τbx and τby are the stress boundary conditions.

Taking all of these expressions together we can write the integrated momentum equations
as,

∂ ∂ 1 h ∂2τ ∂ h ∂2τ ∂b h ∂
(hu) + (hu2 + gh2 − 2
) + (huv − ) = −gh − (PA )
∂t ∂x 2 ρ ∂x ∂y ρ ∂x∂y ∂x ρ ∂x
Z
1 η ∂p 1 ∂2τ ∂2τ
− dz + [ + ],
ρ b ∂x ρ ∂s∂x ∂b∂x
∂ ∂ h ∂2τ ∂ 1 h ∂2τ ∂b h ∂
(hu) + (huv − )+ (hv 2 + gh2 − 2
) = −gh − (PA )
∂t ∂x ρ ∂y∂x ∂y 2 ρ ∂y ∂y ρ ∂y
Z
1 η ∂b 1 ∂2τ ∂2τ
− dz + [ + ],
ρ b ∂y ρ ∂s∂y ∂b∂y

∂ ∂ ∂ 1 1 ∂ ∂2τ ∂ ∂2τ
(hw) + (huw) + (hvw) = − p|z=b + [ (h )+ (h )
∂t ∂x ∂y ρ ρ ∂x ∂z∂x ∂y ∂z∂y
∂η ∂ 2 τ ∂η ∂ 2 τ ∂2τ
−( + − 2 )z=η
∂x ∂z∂x ∂y ∂z∂y ∂z
2 2
∂b ∂ τ ∂b ∂ τ ∂2τ
+( + − 2 )z=b ].
∂x ∂z∂x ∂y ∂z∂y ∂z

Assuming that p is insignificant compared to the other terms and assuming that the lateral
stress averages τxx , τxy and τyy are small simplifies the momentum equations and allows us
to ignore the vertical momentum equation. The pertinent equations are then,

24
 ∂h ∂ ∂
+ (hu) + (hv) = 0, (2.36)
∂t ∂x ∂y
 2 2

∂ ∂ 1 2 ∂ ∂b h ∂ 1 ∂ τ ∂ τ
(hu) + (hu2 + gh ) + (huv) = −gh − (PA ) + + ,
∂t ∂x 2 ∂y ∂x ρ ∂x ρ ∂s∂x ∂b∂x
(2.37)
 
∂ ∂ ∂ 1 ∂b h ∂ 1 ∂2τ ∂2τ
(hu) + (huv) + (hv 2 + gh2 ) = −gh − (PA ) + + .
∂t ∂x ∂y 2 ∂y ρ ∂y ρ ∂s∂x ∂b∂x
(2.38)

By neglecting the hydrostatic pressure and normal shear stresses we arrive at the following
two-dimensional single-phase (SWEs),
¯

∂h ∂ ∂
+ (hu) + (hv) = 0, (2.39)
∂t ∂x ∂y
∂ ∂ 1 ∂ ∂b
(hu) + (hu2 + gh2 ) + (huv) = −gh , (2.40)
∂t ∂x 2 ∂y ∂x
∂ ∂ ∂ 1 ∂b
(hu) + (huv) + (hv 2 + gh2) = −gh . (2.41)
∂t ∂x ∂y 2 ∂y

In above equation we have dropped the bar from the averaged quantities for convenience.
The single phase shallow water equations can be deduced from above equations by ne-
glecting the y-component. The 1-D SPSF system with bottom topography is given below,

 ∂h ∂m

 + = 0,
∂t ∂x   (2.42)

 ∂m ∂ m2 g 2 ∂b
 + + h = −gh ,
∂t ∂x h 2 ∂x

where g is gravitational acceleration constant, flow height is denoted by h , m = hu denotes

the momentum, u is the velocity and b = b(x), x ∈ R is the height of the bottom from a
given level. For steady state solutions, the equilibria of the system satisfy

u = 0 and h + b = constant. (2.43)

25
In compact form, the Eq. (2.42) can be written as

∂Wl ∂Fl
+ = τl , l = 1, 2 , (2.44)
∂t ∂x

where W1 = h, W2 = hu and the fluxes are given as

m2 g 2 (w2 )2 g
F1 = m = w2 , F2 = + h = + (w1 )2 . (2.45)
h 2 w1 2

Moreover,
∂b ∂b
τ1 = 0, τ2 = −gh = −gw1 . (2.46)
∂x ∂x

Hereafter, Wl , l = 1, 2 will be referred as the conserved variables. The quasi-linear form of

the equations is given by

∂W ∂W ∂b
+ A(W) = −gh , (2.47)
∂t ∂x ∂x

where W = (W1 , W2 )T and A(W) is the Jacobian matrix whose element at the k-th row
∂Fk
and l-th column is ∂Wl
for k, l = 1, 2, and the functions Fk (W1 , W2 ) are defined in Eq.
(2.45). The matrix A(W) is given as
 
0 1
A(W) =  w2
, (2.48)
− w22 + gw1 2 ww12
1

system is defined over

Ω = {w ∈ R2 ; h ≥ 0, u ∈ R}. (2.49)

The system has eigenvalues λ± = u ∓ c and the corresponding right eigenvectors r± are

26
given by
 
1
r± =  . (2.50)
u∓c

√
Here, c = gh and the system is strictly hyperbolic for h > 0.

2.4 Derivation of Two-Phase Shallow Flow Model

In this section we have derived and analysed multi-layer SWEs [34]. The multi-layer shallow
water equations are derived from the integration of the same equations as the single-layer
shallow water equations, the inviscid Navier-Stokes equations which consists of equations
for,

∂u
incompressibility∂x
+ ∂w
∂z
= 0,
mass conservation ∂ρ
∂t
∂
+ ∂x (uρ) + ∂z ∂
(wρ),
∂ ∂ ∂ ∂P
x-momentum conservation ∂t (ρu) + ∂x (ρu2 ) + ∂z
(ρuw) = ∂x
], and
∂ ∂ ∂
z-momentum conservation ∂t (ρw) + ∂x (ρuw) + ∂z
(ρw 2 ) = − ∂P
∂z
− ρg,

in one-space dimension. In the multi-layer case the velocity u and density ρ are allowed to
be piece-wise defined through the depth where η defines the interface between each layer.
Integration of the incompressibility condition gives the familiar equation for the depth dy-
namics,

Z η1 Z η1
∂u ∂w ∂ ∂ ∂
( + )dz = udz − (η1 )u|η1 + (η2 )u|η2 + w|η1 − w|η2 ,
η2 ∂x ∂z ∂x η2 ∂x ∂x
∂ ∂
= (hs us ) + (η1 − η2 ),
∂x ∂t
∂ ∂
= (hs ) + (hs us ) = 0.
∂t ∂x

Integration of the mass conservation can be done using the conservative form ρt + (uρ)x +

27
(wρ)z = 0 or through the equation above. In both cases, the integration yields

∂ ∂
(ρs hs ) + (ρs hs us ) = 0,
∂t ∂x

which can be used instead of the continuity equation above. This can be modified such
that the derivatives only exists on ρs as in,

∂ρs ∂ρs
+ us = 0.
∂t ∂x
Integration through the bottom layer gives
Z η2 Z η2
∂u ∂w ∂ ∂ ∂b
( + )dz = udz − (η2 )u|η2 + u|b + w|η2 − w|b ,
b ∂x ∂z ∂x b ∂x ∂x
∂ ∂
= (hf uf ) + (η2 ),
∂x ∂t
∂ ∂
= (hf ) + (hf uf ) = 0,
∂t ∂x
and as before
∂ ∂
(ρf hf ) + (ρf hf uf ) = 0 and
∂t ∂x
∂ ∂
(ρf ) + uf (ρf ) = 0.
∂t ∂x

In order to integrate the horizontal-momentum equations, we first need to find the pressure
as a function of depth. Looking at the vertical-momentun equation, integrating in depth
from the top surface η1 to an intermediate surface coordinate ξ and making the assumption
that the advective terms for w are negligible we have,
Z η1 Z η1
Pz dz = ρgdz. (2.51)
ξ ξ

Integrating the left hand side we have

Z η1
Pz dz = P (x, η1 , t) − P (x, ξ, t) = PA (x, t) − P (x, ξ, t), (2.52)
ξ

28
where PA (x, t) is the surface pressure at η1 . The right hand side of (2.4) gives a different
result depending on where ξ is in relation to the internal surface η2


Z η1  ρ (x, t)g(η − ξ),
1 1 η2 < ξ < η1 ,
ρgdz = (2.53)
ξ  ρ (x, t)gh (x, t) + ρ (x, t)g(η − ξ), b < ξ < η2 .
1 1 2 2

Integrating the horizontal-momentum equation through the top layer leads to,
Z η1 Z η1 Z η1
∂ ∂ 2 ∂ ∂ ∂ ∂ ∂
[ (ρu) + (ρu ) + (ρuw)]dz = ρudz − (η1 )ρu|η1 + (η2 )ρu|η2 + ρu2 dz
η2 ∂t ∂x ∂z ∂t η2 ∂t ∂t ∂x η2
∂ ∂
− (η1 )ρu2 |η1 + (η2 )ρu2 |η2 + ρuw|η1 + ρuw|η2 ,
∂x ∂x
∂ ∂
= (ρs us ) + (ρs hs u2s ),
∂t ∂x

for the left hand side. Integration of the pressure gradient leads to
Z η1 Z η1
∂P ∂ ∂ ∂
− dz = − P dz + (η1 )P |η1 − (η2 )P |η2 ,
η2 ∂x ∂x η2 ∂x ∂x
∂ 1 ∂
= − ( gρs (η1 ) − η2 )2 ) − (hf + b)ρs g(η1 − η2 ),
∂x 2 ∂x
1 ∂
= −( ρs gh2s ) − ρs ghs (hf + b).
2 ∂x

This leads to the integrated x-momentum equation for the top layer i.e

∂ ∂ 1 ∂
(ρs hs us ) + (ρs hs u2s + ρs gh2s ) = −ρs ghs (hf + b). (2.54)
∂t ∂x 2 ∂x

For the second layer, the integration of the left hand side is identical but the pressure
integration is more complicated. Integrating the pressure leads to,
Z η2 Z η2
∂P ∂ ∂ ∂b
− dz = − P dz + (η2 )P |η2 − P |b . (2.55)
b ∂x ∂x b ∂x ∂x

29
Pressure integrated through the depth gives,
Z η2 Z η2
∂ ∂
− P dz = − (gρs hs + gρf (η2 − z))dz,
∂x b ∂x b
∂ 1
= − (gρs hs hf + gρf (η2 − b)2 ),
∂x 2
∂ 1
= − (gρs hs hf + gρf h2f ).
∂x 2

On evaluating the pressure gives,

∂ ∂
(η2 )P |η2 = ρs ghs (hf + b),
∂x ∂x
∂b ∂b
P |b = (ρs ghs + ρf ghf ) .
∂x ∂x

This gives us the following form for the pressure integral through the second layer,
Z η2
∂P ∂ 1 ∂ ∂b
− dz = − (gρs hs hf + gρf h2f ) + ρs ghs (hf + b) − (ρs ghs + ρf ghf ) .
b ∂x ∂x 2 ∂x ∂x
(2.56)

Finally putting this integration with left hand side integration for the bottom layer, we
obtain the following depth integrated conservation law for the bottom layer,

∂ ∂ 1 ∂ ∂b
(ρf hf uf ) + (ρf hf u2f + ρf gh2f ) = −ghf (ρs hs ) − gρf hf . (2.57)
∂t ∂x 2 ∂x ∂x

The new system of equations for the one-dimensional, two-layer, depth-integrated flow with
variable horizontal density is,

∂ ∂
(hs ) + (hs us ) = 0, (for continuity)
∂t ∂x
∂ ∂
(ρs hs ) + (hs us ρs ) = 0, (for mass conservation)
∂t ∂x
∂ ∂ 1 ∂
(ρs hs us ) + (ρs hs u2s + ρs gh2s ) = −ρs ghs (hf + b), (for momentum conservation)
∂t ∂x 2 ∂x

30
in the top layer and,

∂ ∂
(hf ) + (hf uf ) = 0, (for continuity)
∂t ∂x
∂ ∂
(ρf hf ) + (hf uf ρf ) = 0, (for mass conservation)
∂t ∂x
∂ ∂ 1 ∂ ∂b
(ρf hf uf ) + (ρf hf u2f + ρf gh2f ) = −ghf (ρs hs ) − gρf hf , (for momentum conservation)
∂t ∂x 2 ∂x ∂x

in the bottom layer. A common approximation used for the multi-layer SWEs is to assume
that the density is constant inside each of the layer. This gives us system of four equations
in one dimension written as

 ∂ ∂

 (hs ) + (hs us ) = 0,

 ∂t ∂x



 ∂ ∂ 1 ∂
 (hs us ) + (hs u2s + gh2s ) = −ghs (hf + b),
∂t ∂x 2 ∂x (2.58)

 ∂ ∂

 (hf ) + (hf uf ) = 0,

 ∂t ∂x


 ∂ (hf uf ) + ∂ (hf u2 + 1 gh2 ) = −rghf ∂ (hs ) − ghf ∂b ,

f
∂t ∂x 2 f ∂x ∂x

where we have ignored the drag forces, b is the bottom topography, g is the gravitational
constant and

ρf
γ= . (2.59)
ρs

In compact form the system can be written as,

∂Wl ∂Fl
+ = τl , l = 1, 2, 3, 4 , (2.60)
∂t ∂x

where

W1 = hs , W2 = ms = hs us , W3 = hf , W4 = mf = hf uf , (2.61)

31
and the fluxes are

F1 = ms = w2 ,
m2s g 2 g (w2 )2 g g
F2 = + hs + (1 − γ) hs hf = + (w1 )2 + (1 − γ) w1 w3 ,
hs 2 2 w1 2 2
F3 = mf = w4 ,
m2f g (w4 )2 g
F4 = + h2f = + (w3 )2 . (2.62)
hf 2 w3 2

Moreover,

∂hf ∂b ∂w3 ∂b
τ1 = 0 = τ3 , τ2 = −γghs− ghs = −γgw1 − gw1 ,
∂x ∂x ∂x ∂x
∂hs ∂b ∂w1 ∂b
τ4 = −γghf − ghf = −γgw3 − gw3 . (2.63)
∂x ∂x ∂x ∂x

The conservative and non-conservative fluxes of the system are Fl and τl . We can write
momentum of the mixture mm = ms + γmf as,

∂ ∂ ∂b
mm + Fm (W) = −g (hs + γhf ) , (2.64)
∂t ∂x ∂x

where W = (W1 , W2 , W3 , W4 )T and Fm (W) is defined as

Fm (W) =F2 (W) + γF4 (W) + γghs hf

m2s m2f g 1+γ
= +γ + (h2s + γh2f ) + g hs hf . (2.65)
hs hf 2 2

In the above equation, F2 (W) and F4 (W) are the second and fourth components of Fl (W)
defined in Eq. (2.62). We write the system in the form

∂W ∂W
+ A(W) = ψ b (W) , (2.66)
∂t ∂x

32
where, ψ b (W) = (0, −ghs ∂x b, 0, −ghf ∂x b)T , A(W) (Jacobian matrix) is

 
0 1 0 0
 2 
 w4 
− w2 + gw1 + g 1−γ
2
w 3 2 w2
g 1−γ
2
w 1 0 

A(W) =  1 w1 . (2.67)

 0 0 0 1 
 
w2
0 0 − w42 + gw3 2 w
w3
4
3

System is defined over:

Ω = {W ∈ R4 ; hs , hf ≥ 0 and us , uf ∈ R}, (2.68)

hs
or, equivalently, in terms of h = hs + hf and φ = hs +hf
,

Ω = {W ∈ R4 ; h ≥ 0, φ ∈ [0, 1], and us , uf ∈ R}. (2.69)

2.4.1 Steady States

The equilibrium conditions for the system (2.60) are

hs us = hf uf = constant, (2.70)

u2f
+ g (hs + hf + b) = constant, (2.71)
2

 
u2f u2 1−γ
hs ∂x − s +g hf = constant. (2.72)
2 2 2

The above equations becomes (hs + hf + b) = constant, hs /hf = constant, for a particular
hs
case of steady state i.e.(us = uf = 0). As h = hs + hf and φ = hs +hf
, we must have,

h + b = constant and φ = constant (2.73)

33
2.5 The Two-Dimensional Two-Phase Shallow Gran-
ular Flow Model

The two-dimensional TPSF model has the following form

∂ ∂ ∂
(hs ) + (mxs ) + (my ) = 0 , (2.74)
∂t ∂x ∂y s
 
∂ x ∂ (mxs )2 1 2 1−γ ∂ mxs mys
(ms ) + + ghs + g hs hf + ( ) = 0, (2.75)
∂t ∂x hs 2 2 ∂y hs
 
∂ y ∂ mxs mys ∂ (mys )2 1 2 1−γ
(m ) + ( )+ + ghs + g hs hf = 0 , (2.76)
∂t s ∂x hs ∂y hs 2 2
∂ ∂ ∂
(hf ) + (mxf ) + (myf ) = 0 , (2.77)
∂t ∂x ∂y
 x 2  y
∂ x ∂ (mf ) 1 2 ∂ mxf mf
(m ) + + ghf + ( ) = 0, (2.78)
∂t f ∂x hf 2 ∂y hf
!
x y y 2
∂ ∂ mf m f ∂ (m f ) 1
(my ) + ( )+ + gh2f = 0 . (2.79)
∂t f ∂x hf ∂y hf 2

Here, hs and hf are the heights of solid and fluid phases respectively. Also, mxs = hs us and
mys = hs vs are respectively the momentums of solid phase in x and y directions and us and
vs are the components of fluid velocity vector u. Similarly, mxf = hf uf and myf = hf vf
are respectively the momentums of fluid phase in x and y directions and uf and vf are the
components of fluid velocity vector u.The gravitational constant is represented by g and
ρf
γ= ρs
. In compact form the above system can be written as,

∂Wl ∂Fl
+ = 0, l = 1, 2, · · · , 6 (2.80)
∂t ∂x

where

W1 = hs , W2 = mxs = hs us , W3 = mys = hs vs , (2.81)

W4 = hf , W5 = mxf = hf uf , W6 = myf = hf vf , (2.82)

34
fluxes are

F1 = mxs = w2 ,
(mxs )2 g 2 g (w2 )2 g g
F2 = + hs + (1 − γ) hs hf = + (w1 )2 + (1 − γ) w1 w4
hs 2 2 w1 2 2
x y
m m w2 w3
F3 = s s = , and F4 = mxf = w5 ,
hs w1
x 2
(mf ) g 2 (w5 ) 2
g 2
mxf myf w5 w6
F5 = + hf = + (w4 ) F6 = = (2.83)
hf 2 w3 2 hf w4

and

mxs mys w2 w3
g1 = mys = w3 , g2 = = ,
hs w1
(mys )2 g 2 g (w3 )2 g g
g3 = + hs + (1 − γ) hs hf = + (w1 )2 + (1 − γ) w1 w4 ,
hs 2 2 w1 2 2
x y
mf mf w5 w6
g4 = myf = w6 , g5 = = ,
hf w4
(myf )2 g 2 (w6 )2 g
g6 = + hf = + (w4 )2 . (2.84)
hf 2 w4 2

35
Chapter 3

The CESE Scheme

This chapter introduces the derivation of CESE scheme for SPSF and TPSF models.

3.1 Numerical Scheme

CESE scheme is the most famous and frequently used scheme among the present schemes
for solving system of hyperbolic conservation laws. Since its inception, a lot of progress
has been made in development and applications of CESE scheme over the past two decades
[23, 24, 25, 27]

3.2 The CESE Method for Single-Phase Shallow Flows

The scheme was made to solve the conservation laws which governs the motion of fluid.
It is based on a different concept, differs from the familiar methods e.g. the FDM, FVM
and FEM. CESE method is made to overcome the disadvantages of these methods [23, 24,
25, 27]. In this section, the one-dimensional CESE method of Chang [23] is extended for
single-phase shallow flow model (2.42). In component form, the Eq. (2.42) can be rewritten
as
∂Wl ∂Fl
+ = 0, l = 1, 2 . (3.1)
∂t ∂x

36
Let x0 = t,x1 = x be the coordinates of 2D space (E2 ) and hl = (Wl , Fl )T , l = 1, 2.

By applying Gauss divergence theorem, Eq.(3.1) becomes

I
hl · ds = 0, l = 1, 2 , (3.2)
S(V )

where (a) hl = (Wl , Fl ), l = 1, 2, i.e. for each l = 1, 2, the components of the vector hl in
the t- and x-directions are Wl and Fl , respectively and (b) ds = dσn. Eq. (3.2) is applied
at space-time domain, known as CE (conservation element). Integration is being done by
the use of SEs (solution elements). When we talk about SE, we suppose that flow variables
are smooth. Thus, variables of flow can be described having prescribed order of accuracy.

We denote by Ω as the set of mesh points (j, n) in E2 space where n = 0, 21 , 1, 32 , · · · and

for each n, j = 0, ± 12 , ±1, ± 32 , · · · , for each mesh point (j, n), we have a solution element
associated with it shown in the Fig. 3.1 which is denoted by dash curve. It does not matter
what is the exact size of the neighborhood. Now Wl (x, t), Fl (x, t) and hl (x, t) are approx-
imated by Wl∗ (x, t; j, n), Fl∗ (x, t; j, n) and h∗l (x, t; j, n) respectively as follows.

aWl∗ (x, t; j, n) = (Wl )nj + (Wlt )nj (t − tn ) + (Wlx )nj (x − xj ) , (3.3)

and

Fl∗ (x, t; j, n) = (Fl )nj + (Flt )nj (t − tn ) + (Flx )nj (x − xj ). (3.4)

By employing chain rule, we obtain

4
X
(Flx )nj = (Fl,m )nj (Wmx )nj , (3.5)
m=1

37
 4
X
(Flt )nj = (Fl,m )nj (Wmt )nj , (3.6)
m=1

where

∂Fl
Fl,m = , l, m = 1, 2 . (3.7)
∂Wm

Jacobian matrix formed by Fl,m l, m = 1, 2 is given by Eq. (2.48). Note that (Wl )nj ,(Wl x)nj
∂Wl
and (Wl t)nj are constant in SE(j, n). These are numerical analogues of values of Wl , ∂x
∂Wl
and ∂t
at (xj , tn ) respectively. Since hl = (Wl , Fl )T , we can write

h∗l (x, t : j, n) = (Wl∗ (x, t : j, n), Fl∗ (x, t : j, n)) . (3.8)

Moreover, due to the smoothness assumption of variables, Wl∗ (x, t; j, n) and Fl∗ (x, t; j, n)
satisfy the Eq.(23), i.e.,

∂Wl∗ (x, t; j, n) ∂Fl∗ (x, t; j, n)

+ = 0. (3.9)
∂t ∂x

Using Eq. (3.3) and (3.4), Eq. (3.9) is equivalent to

(Wlt )nj = −(Flx )nj . (3.10)

We notice that (Fl )nj are functions of (Wl )nj , (Flx )nj are functions of (Wl )nj and (Wlx )nj , (Flt )nj
are functions of (Wl )nj and (Wlt )nj .
Consider the Fig. 3.1,

38
 (a)
1 1
j− 3 j − 1j − 2
j j+ 2 j+1 j+ 3
2 2
n+1
1
n+ 2
∆t/2
n
∆t/2
1
n− 2

t n−1
x ∆x/2 ∆x/2
(j,n) (j,n) (j,n)
B A A F B F
A
C D D E C E
D
(j−1/2,n−1/2) (j+1/2,n−1/2) (j−1/2,n−1/2) (j+1/2,n−1/2)
(b) (c)

Figure 3.1: Space-time staggered grid near SE (j, n)

The rectangular non-overlapped areas are CEs. The conservation element at (j, n) is shown
as the dotted rectangles. The integral representation of Equation. (3.1) is given as

Z
Wl∗ dx − Fl∗ dt = 0 , l = 1, 2 (3.11)
∂Ω

The above equation represents the global flux conservation for whole domain Ω. Due to
local flux conservation property, this equation is still satisfied for any conservation element.
Integrating at mesh points (j ± 12 , n − 12 ) one can obtain

n− 1 ∆x h n− 1
i
(Wl )nj =(Wl )j± 12 ∓ (Wlx )j± 12 + (Wlx )nj (3.12)
2 4 2

∆t h 1 i ∆t2 h i
n− 2 n n− 12 n
∓ (Fl )j± 1 − (Fl )j ∓ (Flt )j± 1 + (Flt )j .
∆x 2 4∆x 2

39
The summation gives

1h n− 1 n− 1 n− 1 n− 1
i
(Wl )nj = (Wl )j− 12 + (Wl )j+ 12 + (Sl )j− 12 − (Sl )j+ 12 , (3.13)
2 2 2 2 2

where
n− 1 ∆x n− 1 ∆t n− 1 ∆t2 n− 1
(Sl )j± 12 = (Wlx )j± 12 + (Fl )j± 12 + (Flt )j± 12 . (3.14)
2 4 2 ∆x 2 4∆x 2

Note that, T1 = 0 = T3 = 0. Moreover,

n− 21 n− 12 ∆t  n− 21 n− 21
h
n− 12 n− 21
i
(T2 )j+ 1 − (T2 )j− 1 = γg (hs )j− 1 + (hs )j+ 1 (hf )j+ 1 − (hf )j− 1 , (3.15)
2 2 2∆x 2 2 2 2

and
n− 12 n− 12 ∆t  n− 12 n− 21
h
n− 21 n− 21
i
(T3 )j+ 1 − (T3 )j− 1 =g (hf )j− 1 + (hf )j+ 1 (hs )j+ 1 − (hs )j− 1 . (3.16)
2 2 2∆x 2 2 2 2

For slopes of conserved variables, we can use limiting formulations and thus can suppressed
the numerical oscillations near a discontinuity.

(Wlx )nj = Ul ((Wlx− )nj , (Wlx+ )nj ; α) , l = 1, 2 . (3.17)

Now we have, α ≥ 0, a constant (normally α = 1 or α = 2) and

|x+ |α x− + |x− |α x+
Ul (x− , x+ ; α) = . (3.18)
|x+ |α + |x− |α

Moreover,
(Wl′)nj+ 1 − (Wl )nj (Wl )nj − (Wl′ )nj− 1
(Wlx+ )nj = 2
, (Wlx− )nj = 2
l = 1, 2 (3.19)
∆x/2 ∆x/2

40
and
n− 1 ∆t n− 1
(Wl′)nj± 1 = (Wl )j± 12 + (Wlt )j± 12 l = 1, 2 . (3.20)
2 2 2 2

Eqs. (3.13) and (3.20) constitute the CESE solver for SPSF model. The steep wave fronts
are treated by approximating (Wlx )nj through averaging procedure described in above. This
completes the derivation of proposed solver.

3.3 The CESE Method for Two-Phase Shallow Flows

The CESE Scheme was made to solve conservation laws which governs the motion of
fluid. It is based on a different concept, differs from the familiar methods e.g. FDM,
FVM and FEM. CESE method is made to overcome the disadvantages of these methods
[23, 24, 25, 27]. In this section, the one-dimensional CE/SE method of Chang [23] is ex-
tended for two-phase shallow flow model (2.60). Coordinates for E2 are the same as in the
case for SPSF.

By applying Gauss divergence theorem, we can write Eq.(2.60) as,

I Z
hl · dS = τl dV, l = 1, 2, 3, 4 , (3.21)
S(V ) V

where (a) hl = (Wl , Fl )T , l = 1, 2, 3, 4, and (b) dS = dσn where dσ is area, n is unit out-
ward normal vector. Eq. (3.21) is applied at space-time domain, known as conservation
element (CE). Integration is being done by the use of SEs (solution elements). When we
talk about SE, we suppose that flow variables are smooth. Thus variables of flow can be
described having prescribed order of accuracy.

We denote by Ω as the set of mesh points (j, n) in E2 space where n = 0, 12 , 1, 32 , · · · and

for each n, j = 0, ± 12 , ±1, ± 32 , · · · . For each mesh point (j, n), we have a solution element

41
(SE) associated with it shown in the Figure 3.1, denoted by a dash curve. It does not
matter what is the exact size of the neighborhood. Now Wl (x, t), Fl (x, t) and hl (x, t) are
approximated by Wl∗ (x, t; j, n), Fl∗ (x, t; j, n) and h∗l (x, t; j, n) respectively as follows.

Wl∗ (x, t; j, n) = (Wl )nj + (Wlt )nj (t − tn ) + (Wlx )nj (x − xj ) , (3.22)

and

Fl∗ (x, t; j, n) = (Fl )nj + (Flt )nj (t − tn ) + (Flx )nj (x − xj ) . (3.23)

Moreover,
τl∗ (x, t; j, n) = τ (Wl∗ , Wlx∗ ) , (3.24)

By employing chain rule, we obtain

4
X
(Flx )nj = (Fl,m )nj (Wmx )nj , (3.25)
m=1

4
X
(Flt )nj = (Fl,m )nj (Wmt )nj (3.26)
m=1

where
∂Fl
Fl,m = , l, m = 1, 2, 3, 4 . (3.27)
∂Wm

Jacobian matrix formed by Fl,m , l, m = 1, 2, 3, 4 , is given by Eq. (2.67). Note that (Wl )nj ,
(Wlx )nj and (Wlt )nj are constant in SE(j, n). Since hl = (Wl , Fl )T , we write

h∗l (x, t : j, n) = (Wl∗ (x, t : j, n), Fl∗ (x, t : j, n))T . (3.28)

Moreover, due to the smoothness assumption of variables, Wl∗ (x, t; j, n) ,Fl∗ (x, t; j, n) and

42
τl∗ (x, t : j, n) satisfy the Eq. (2.60), i.e.,

∂Wl∗ (x, t; j, n) ∂Fl∗ (x, t; j, n)

+ = τl∗ (x, t; j, n) . (3.29)
∂t ∂x

Using Eqs. (3.22) and (3.23), Eq. (3.29) is equivalent to

(Wlt )nj = −(Flx )nj + τl ((Wl )nj , (Wlx )nj ) . (3.30)

We notice that (Fl )nj are functions of (Wl )nj , (Flx )nj are functions of (Wl )nj and (Wlx )nj , (Flt )nj
are functions of (Wl )nj and (Wlt )nj , and τl are functions of (Wl )nj and (Wlx )nj .
Let E2 be divided into non-overlapping rectangular regions (see Figure 3.1(a)) referred to
as conservation elements(CEs).
Therefore we obtain,

I Z
∗
h l · dS = τl dV , l = 1, 2, 3, 4 . (3.31)
S(CE± (j,n)) CE± (j,n)

Clearly

I Z
∗
h l · dS = τl dV , l = 1, 2, 3, 4 , (j, n) ∈ Ω , (3.32)
S(CE(j,n)) CE(j,n)

comes from Eq. (3.31).

By using Eq. (3.31) along with Eqs. (3.22), (3.23), (3.24) and (3.28), one obtains

n− 1 ∆x h n− 1
i
(Wl )nj =(Wl )j± 12 ∓ (Wlx )j± 12 + (Wlx )nj (3.33)
2 4 2

∆t h n− 12 n
i ∆t2 h
n− 12 n
i h
n− 12 n
i
∓ (Fl )j± 1 − (Fl )j ∓ (Flt )j± 1 + (Flt )j ± (τl )j± 1 − (τl )j .
∆x 2 4∆x 2 2

43
The summation gives

1h n− 1 n− 1 n− 1 n− 1 n− 1 n− 1
i
(Wl )nj = (Wl )j− 12 + (Wl )j+ 12 + (Sl )j− 12 − (Sl )j+ 12 + (τl )j+ 12 − (τl )j− 12 , (3.34)
2 2 2 2 2 2 2

where
n− 1 ∆x n− 1 ∆t n− 1 ∆t2 n− 1
(Sl )j± 12 = (Wlx )j± 12 + (Fl )j± 12 + (Flt )j± 12 . (3.35)
2 4 2 ∆x 2 4∆x 2

Note that, τ1 = 0 = τ3 = 0. Moreover,

n− 1 n− 1 ∆t  n− 1 n− 1
h
n− 1 n− 1
i
(τ2 )j+ 12 − (τ2 )j− 12 =γg (hs )j− 12 + (hs )j+ 12 (hf )j+ 12 − (hf )j− 12
2 2 2∆x 2 2 2 2

∆t  n− 12 n− 12
h
n− 12 n− 21
i
+g (hs )j− 1 + (hs )j+ 1 (b)j+ 1 − (b)j− 1 , (3.36)
2∆x 2 2 2 2

and
n− 21 n− 12 ∆t  n− 12 n− 21
h
n− 12 n− 21
i
(τ3 )j+ 1 − (τ3 )j− 1 =g (hf )j− 1 + (hf )j+ 1 (hs )j+ 1 − (hs )j− 1
2 2 2∆x 2 2 2 2

∆t  h i
n− 12 n− 12 n− 21 n− 12
+g (hf )j− 1 + (hf )j+ 1 (b)j+ 1 − (b)j− 1 . (3.37)
2∆x 2 2 2 2

To achieve the steady state conditions described in Eq. (2.71), we have discretized the
bottom topography term in the same manner as given in [17]. Whenever Riemann data
corresponds to equilibrium, the steady state condition holds i.e. when (h + b)ni = (h + b)ni+1 ,
φni = φni+1 , (us )ni = (us )ni+1 = (uf )ni = (uf )ni+1 = 0.

For slopes of conserved variables, we can use limiting formulations and thus can suppressed
the numerical oscillations near a discontinuity.

(Wlx )nj = Ul ((Wlx− )nj , (Wlx+ )nj ; α) , l = 1, 2, 3, 4 . (3.38)

44
Now we have, α ≥ 0, a constant (normally α = 1 or α = 2) and

|x+ |α x− + |x− |α x+
Ul (x− , x+ ; α) = . (3.39)
|x+ |α + |x− |α

Moreover,

(Wl′ )nj+ 1 − (Wl )nj (Wl )nj − (Wl′ )nj− 1

(Wlx+ )nj = 2
, (Wlx− )nj = 2
, (3.40)
∆x/2 ∆x/2

and

n− 1 ∆t n− 1
(Wl′ )nj± 1 = (Wl )j± 12 + (Wlt )j± 12 , l = 1, 2, 3, 4 . (3.41)
2 2 2 2

Eqs. (3.34) and (3.38)-(3.41) constitute the CESE solver for TPSF model. The steep wave
fronts are treated by approximating (Wlx )nj through averaging procedure described above.
This completes the derivation of proposed solver.
The method can be easily extended to solve two-dimensional shallow flow models by using
the approach of Zhang et al. [36] on regular rectangular grids. In this manuscript, we omit
further details about the two-dimensional CESE method. For more details, the reader is
referred to the articles of Zhang et al. [36] and Qamar et al. [37].

45
Chapter 4

Case Studies

Here, some case studies are carried out for the considered models to validate our numerical
scheme. For comparison the models were solved by staggered central scheme [33]. In all
1-D test problems, the reference solutions were obtained on 1000 grid points whereas the
reference solution in the 2-D problems were obtained on 500 × 500 grid points.

4.1 Single-Phase Shallow Flows

4.1.1 Test Problem 1: Rarefaction into Vacuum

Consider a test problem already discussed in [20, 44], which depicts initial dry bed state.
Riemann data for right and left states is hr = 1, hl = 0 (hr = 0, because of the dry right
bed). Also ur = 0 = ul and the constant of gravitation (g) is 1. At time t = 1 momentum
m = hu and the flow height h are shown in Fig. 4.2. An accurate agreement is observed
by comparing the CESE and central schemes.

4.1.2 Test Problem 2: Formation of Dry bed

This test problem was suggested in [20] by Toro. Generation of a dry bed region is included
in it. The Riemann data for the problem is hr = hl = 0.1 and ul = −3, ur = 3 and the
gravitational constant is g = 9.81. Numerical results for height h and momentum m = hu

46
by CESE as well as central schemes are displayed in Fig. 4.2. It is observed that the flow
height has some negative values, although, the computation is carried out for a very small
CFL number.

4.1.3 Test Problem 3: Test with Variable Bottom Topography

This problem was also considered in [45]. Initially, the total height is h = 1 ,velocity u = 0
everywhere and b(x) 6= 0. The bottom topography is taken as,

2 2
b(x) = 0.2e−0.5(x+1) + 0.3e−(x−1.5) . (4.1)

The computational domain is [-10,10]. The steady solutions obtained by using CESE
scheme. The results are displayed in Fig.4.3. It was observed that CESE gives the accurate
approximation to the exact solution u = 0.

4.2 Test Problems of Two-Dimensional Single-Phase

Shallow Flows
Results of 2D SPSF model are shown here. Numerical results of the two schemes were
compared with each other. We discretise the computational domain [0, 2] × [0, 2] into
100 × 100 grid points. The reference solutions are computed on 500 × 500 grid points.

4.2.1 Test Problem 4: Two Separated Fluids

This problem was considered by Qamar et al. in [46] for shallow water magnetohydrody-
namic model. Here, we consider the same data in which magnetic effects are neglected.
The domain is given by [x, y] ∈ [0, 2] × [0, 2] which is a square region with extrapolation
boundaries. The square region contains a circular region of radius r = 0.3 in its center.
The data for both regions are

h = 10 u = (0, 0), if r ≤ 0.3 ,

(4.2)
h = 1.0 u = (0, 0) if r > 0.3 .

47
 Height at t=1

1 CE/SE
central
0.8 reference

0.6
h

0.4

0.2

0
−4 −2 0 2 4
x−axis
Momentum at t=1
0.3
CE/SE
0.25 central
reference
0.2
hu

0.15

0.1

0.05

0
−4 −2 0 2 4
x−axis

Figure 4.1: Rarefaction into vacuum for SPSF.

48
 Height at t=1

0.1 central
CE/SE
0.08 reference

0.06
h

0.04

0.02

−5 0 5
x−axis
Momentum at t=1
0.3

0.2

0.1
hu

−0.1
central
−0.2 CE/SE
reference
−0.3
−5 0 5
x−axis

Figure 4.2: Dry bed formation for SPSF at t=1.

49
 Results with Central scheme at t=12

1
h+b
0.8 u

0.6
y−axis

0.4

0.2

−0.2
−6 −4 −2 0 2 4 6
x−axis
Results with CESE scheme at t=12

1
h+b
0.8 u

0.6
y−axis

0.4

0.2

−0.2
−6 −4 −2 0 2 4 6
x−axis

Figure 4.3: Problem 3: Results of variable topography at t=8 and t=12.

50
The gravitational constant is g = 1. Simulation results are shown in Fig. 4.4 at t = 0.15.
In that figure the cross-sectional solution of CESE scheme along y = 1 are compared with
the central scheme and reference solutions. Slowly moving Alfven waves which are located
between rarefaction and shock waves are visible in the tangential component of velocity.
Both schemes have comparable results.

4.2.2 Test Problem 5: Rotor Like Problem

This problem is also considered in [46], but here we look for solution in which magnetic
effects are neglected. The domain for the problem is a square with a circular section with
radius r = 0.2 centered in (1, 1). Also the boundaries are extrapolated. The initial data
for the problem are

h = 10, u = (−y, x), if r ≤ 0.2 ,

(4.3)
h = 1, u = (0, 0), if r > 0.2 ,

1
where, r=(x2 +y 2) 2 and the gravitational constant g=1. The final simulation time t is 0.15.
Results are shown in Fig. 4.5. Both CESE and central schemes are in close agreement.
However, CESE scheme solution is more closer to the reference solution.

4.3 Two-Phase Shallow Flows

Here, we discuss numerical simulations of the TPSF model (2.58). We set γ = 1/2 for all
the test problems. The CESE scheme is applied on all the six test problems. The results
are compared with central scheme.

No Dry Bed Zones Tests

Results of these problems were considered first. Drag forces are absent in these problems.

51
 Central
3
CESE
2 reference
2
h

1.5
1

h
1
0
2
2
1 0.5
1
y−axis 0 0 x−axis 0.5 1 1.5
x−axis

1.5
Central
2
1 CESE
1 reference
y−velocity

0.5
0
y−velocity

0
−1
−0.5
−2
2
2 −1
1
1
−1.5
y−axis 0 0 x−axis 0.5 1 1.5 2
x−axis

Figure 4.4: Problem 4: Two separated fluids results on 100 × 100 grid points at t=0.15.

52
 0
Central
2 CESE
−0.5 reference
0
x−velocity

x−velocity
−1
−2
−1.5
−4
2
2 −2
1
1
y−axis 0 0 x−axis 0.5 1 1.5 2
x−axis

3 Central
4 CESE
2 reference
2
y−velocity

y−velocity

1
0
0
−2
2
2 −1
1
1
y−axis 0 0 x−axis 0.5 1 1.5 2
x−axis

Figure 4.5: Problem 5: Rotor like problem results on 100 × 100 grid points t = 0.15.

53
4.3.1 Test Problem 6: Simple Riemann Problem
Here, the problem is considered that was suggested in [17]. At x = 0 we have initial
discontinuity location. Riemann data for,

(h, φ, us , uf )l = (1.0, 0.8, 0.0, 0.0) (4.4)

(h, φ, us , uf )r = (2.0, 0.4, −0.9, −0.1), (4.5)

We have considered g = 9.81. The results are performed on 200 mesh cells and the
computational domain is [-5,5]. Time t = 0.5, the results are shown in Fig.4.7. We display
results for solid volume fraction φ, height variables h,hs ,hf ,phase velocities us and uf in sub
figures (a)-(c) of Fig 4.7. one-rarefaction, two-shock, and a three-rarefaction was for the
solution of this Riemann problem.The discontinuities poses us of the possible convergence
difficulties. Here we have compared our results with central schemes and a close agreement
was observed. Also our results gives close agreement with the results of [17].

4.3.2 Error Analysis

In order to verify the performance and accuracy of our scheme, we again consider the test
problem in Sect. 4.1. We compute the errors of height of flow h and solid volume fraction φ
by both CESE and central schemes. L1 errors for h and φ on different mesh sizes are given
in tables 4.1 and 4.2 respectively. The comparison shows better performance of CESE
scheme with less error.

4.3.3 Test Problem 7: Rarefaction of the Fluid Constituent

into Vacuum
Here, consider flow (for h > 0) at spatial domain and ∀ t, but the fluid phase (hf = 0) is
characterized by a vacuum zone. The data for the Riemann problem is:

left state data = (1.0, 0.8, 0.0, 0.0) , if x ≤ 0.0 (4.6)

right state data = (1.0, 1.0, 0.0, 0.0) , if x > 0.0 . (4.7)

54
 Error

CESE
0.03 Central

0.025
velocity

0.02

0.015

0.01

0.005

500 1000 1500

N
Error

CESE
0.03
Central
0.025
velocity

0.02
0.015
0.01
0.005

500 1000 1500

N

Figure 4.6: : Error Analysis for height h and solid volume fraction φ.

55
 Table 4.1: L1 -errors for phi Table 4.2: L1 -error velocity
N L1 -error (Central) L1 -error (CESE) N L1 -error (Central) L1 -error (CESE)

50 0.0268 0.3020 50 0.0011 0.0104

100 0.0559 0.016 100 0.0084 0.0062
200 0.0283 0.0081 200 0.0045 0.0032
400 0.0147 0.0041 400 0.0023 0.0016
800 0.0073 0.0016 800 0.0011 0.0007
1600 0.0038 0.0005 1600 0.0004 0.0002

There is solid only on the right. At x = 0, discontinuity is located and we take g = 9.81.
The computation is performed on domain [−5, 5] with 100 mesh cells. Results are shown
in Fig. 4.8. The solution contains 2-rarefaction (across which hf = h(1 − φ) vanishes),
1-shock, and 4-rarefaction occurring in solid material. Figures depicts excellent agreement
of both CESE and central schemes.
Dry bed Zones Tests
Here, we consider test problems which involves dry bed zones. Height h of flow is taken
zero for h < ǫ = 10−5 (tolerance)

4.3.4 Test Problem 8: Spreading of a Granular Mass

In this case, initially uf = us = 0, and the data for height of flow, solid volume fraction is
given as:

 1.0 if x ∈ [−1, 1], 2
h(x, 0) = ϕ(x, 0) = 0.3 + 0.4e−x . (4.8)
 0.0 otherwise,

We set g = 1. Cases of infinitely large inter -phase drag forces and no drag forces, the solu-
tion is computed on 200 mesh cells with computational domain [-10, 10]. In the problem,
phase velocity differences approaches 0, for no drag forces, as h tends to 0. The results by
CESE scheme and comparison with central scheme are approximated and displayed in the
Fig. 4.9. We observe that the hyperbolic conditions are maintained as flow height vanishes.

56
Also, we observe that height dynamics are not effected by drag forces. The behavior of
solid volume fraction φ varies also noticeably.

4.3.5 Test Problem 9: Generation of Dry Bed

Formation of Dry Bed Zone test problems are considered here similar to Riemann prob-
lems already presented for single-phase model. The solution is consisting of 2-opposite
rarefactions and a dry bed zone is generated between them. Simple Riemann data for the
set of test problems is given as:

Test Problems lef tstatedata rightstatedata

i (0.1, 0.4, −3.0, −3.0) (0.1, 0.7, 3.0, 3.0)

ii (0.1, 0.4, 0.0, 0.0) (0.7, 0.7, 6.0, 6.0)
iii (0.2, 0.4, −3.0, −3.0) (0.1, 0.8, 3.0, 3.0)

In all problems, at x = 0 we have the initial interface. Results are performed over 200 mesh
cells and the computational domain is [-5,5]. In all the three problems, infinitely large drag
forces are applied so that the phase velocity difference approaches 0. This condition assures
our solution in the region of hyperbolicity over temporal and spatial domain.
The test problem (i) is similar to problem in Fig. 4.2, for SPSF model except that here
the solid volume fraction φ has an initial discontinuity. At time t = 1 results are reported
in Fig. 4.10. The hight variables h, hs , hf , momenta mf , ms , mm , solid volume fraction φ,
and velocities us , uf are displayed in sub-figures (a − d) of Fig. 4.10. By comparison, we
observe a a better and qualitative agreement between both schemes.
Test problem (ii) is analogous to test problem 1, except that here, we consider an initial
velocities translation as a result left rarefaction is transonic.
In the last problem we have initial discontinuity in the height of flow h and in the solid
volume fraction φ Fig.4.11 and Fig.4.12, show the results for test problem (ii) and test
problem (iii) respectively.

57
4.3.6 Test Problem 10: Interface Propagation
The aim of this problem is to validate the performance of our suggested scheme (CESE) to
capture and propagate the discontinuities over time. This example has been widely used
to validate schemes for multilayer flows [33, 40]. The data for the problem is given as

(hs , hf , us , uf )l = (0.5, 0.5, 2.5, 2.5) , if x ≤ 0.5 , (4.9)

(hs , hf , us , uf )r = (0.55, 0.45, 2.5, 2.5) , if x > 0.5 . (4.10)

ρs
Initially at x = 0.5 discontinuity is located, g = 9.81, and γ = ρf
= 0.98. The computa-
tional domain is [0, 1] which is subdivided into 400 mesh cells for the numerical solution
and 3000 cells for the reference solution of CESE. Fig. 4.13 shows the evolution of interface
between the two layers, free surface displacement and the velocity of each layer. A close
agreement can be observed between the two schemes. However, CESE scheme resolves
better the shock discontinuities. The results also confirms the ability of both schemes to
locate and resolve discontinuities in the solution.

4.4 Test Problem with Variable Bottom Topography

We consider following cases here to check the effectiveness of the CESE scheme.

4.4.1 Perturbation of a Steady State at Rest

This test problem is the extension of Le Veque’s test [39] for SPSF equations having bottom
topography. This problem was also presented by M. Pelanti et. al in[17]. We define Bottom
topography as

 0.25(cos(π(x − 0.5)/0.1) + 1), if |x − 0.5| < 0.1,
b(x) = (4.11)
 0 otherwise.

where small perturbation behavior of steady-state conditions are being observed at rest
(2.73). In the begining, we took a small perturbation of φ, flow depth h:

h(x, 0) = ho + h̃ φ(x, 0) = φo − φ̃ for − 0.6 < x < −0.5, (4.12)

58
with ho = 1, φo = 0.6, and h̃ = φ̃ = 10−3. Fig.4.14 gives the type of initial data for
solid volume fraction φ and total height h + b. We took g = 1 and γ = 1/2, free flowing
boundary conditions were used and computational domain was [−0.9, 1.1]. The solution is
computed at 200 grid points and the reference solution is computed on 3000 mesh cells.

4.4.2 Two-Dimensional Two-Phase Problem

To check the performance of CESE scheme for two-dimensional TPSF model, we consider
a square domain with a circular section with radius r = 0.1 centered in (0, 0). Also the
outflow boundary conditions are considered here. The initial data for the problem are

h = 1 (us , vs ) = (1, 1) (uf , vf ) = (1, 1) φ = 0.8 if r ≤ 0.1 , (4.13)

h = 2 (us , vs ) = (1, 1) (uf , vf ) = (1, 1) φ = 0.4 if r > 0.1 , (4.14)

1
where, r=(x2 + y 2) 2 and the gravitational constant g=1. The final simulation time t is 0.5.
Results are shown in Fig. 4.15.

59
 Height at t=0.5
3 central
h CE/SE
reference
2.5
h
s
h 2

1.5

hf
1

−5 0 5
(a)
Solid Volume fraction φ at t=0.5
0.7
central
0.65 CE/SE
reference
0.6

0.55
φ

0.5

0.45

0.4
−5 0 5
(b)
velocities at t=0.5
1.5

0.5 us
u

−0.5
central
−1 CE/SE
uf reference
−1.5
−5 0 5
(c)

Figure 4.7: Test Problem 5: Comparison of results at t = 0.5 over 200 cells.
60
 Height at t=1
h
1
h
s
0.8 central
CE/SE
h 0.6 reference

0.4
h
f
0.2

0
−5 0 5
(a)
Solid volume fraction φ at t=1

1.1 central
CE/SE
reference
1
φ

0.9

0.8

0.7
−5 0 5
(b)
velocities at t=1

1 central
u
s CE/SE
0.8
reference
u
0.6 f

0.4
u

0.2

−0.2

−4 −2 0 2 4
(c)

Figure 4.8: Test Problem 7: Comparison of results over 100 cells.

61
 Height

0.2

h
0.15
hs
h

0.1 hf

0.05 central
CE/SE
reference
0
−10 −5 0 5 10
(a)
veolicities

1.5 central
CE/SE
1 reference
0.5
0
u

−0.5
u
s
−1
u
f
−1.5

−10 −5 0 5 10
(c)

Figure 4.9: Test Problem 8: Comparison of results for granular mass spreading.

62
 Height at t=1

0.1 central
h CE/SE
0.08 reference

0.06
hs
h

0.04
hf

0.02

−5 0 5
(a)
Momenta at t=1

central
0.2
CE/SE
reference
0.1
hu

0
mm
−0.1
m
s
−0.2
m
f

−5 0 5
(b)
velocities at t=1
3
central
CE/SE
2 reference

1 u
s

0
u

u
f

−1

−2

−3
−5 0 5
(d)

Figure 4.10: Test Problem 9 (i): Comparison of results for generation of dry bed problem.
63
 Height at t=0.5

0.1
h

0.08
hs
h 0.06
hf
0.04
central
0.02 CE/SE
reference

−5 0 5
(a)
Momenta at t=0.5
mm
central
0.5
CE/SE
ms
reference
0.4

0.3
hu

mf
0.2

0.1

0
−5 0 5
(b)
velocities at t=0.5
6
Central
CE/SE
5
reference

4 us

3
u

uf

0
−5 0 5
(d)

Figure 4.11: Test Problem 9 (ii), Comparison of results over domain [-5,5] at time t=1
64
 Height at t=1
0.2 Central
CE/SE
h reference
0.15

h hs

0.1 hf

0.05

−5 0 5
(a)
Momenta at t=1

central
0.2
CE/SE
0.1 reference

0
hu

−0.1 mm

−0.2
ms
−0.3
−0.4
m
f
−0.5
−5 0 5
(b)
velocities at t=1
3
central
2 CE/SE
reference
1

0
u

us
−1
uf
−2

−3
−5 0 5
(d)

Figure 4.12: Test Problem 9 (iii), Results

65 of for height, momenta and velocities.
 Interface evolution Free surface displacement

0.55 CE/SE CE/SE

Central 1.0002 Central
reference reference
0.54
1.0001
0.53
h1

h
1
0.52
0.9999
0.51
0.9998
0.5
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x−axis x−axis
Velocity of lower layer Velocitiy of upper layer
2.5
CE/SE 2.512 CE/SE
2.498 Central Central
2.51
reference reference
2.496
2.508
2.494 2.506
u1

u2

2.492 2.504

2.49 2.502

2.488 2.5
0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1
x−axis x−axis

Figure 4.13: Problem 10 : Results for Interface evolution, total height and velocities.

66
 h+b at t=0.25 φ at t=0.25
1.0005 0.6004
CE/SE CE/SE
reference 0.6002 reference

0.6

0.5998
h+b

φ
0.5996

0.5994

0.5992

0.9995 0.599
−0.5 0 0.5 1 −0.5 0 0.5 1
x−axis x−axis
h+b at t=1.25 φ at t=1.25
1.0005 0.6004
CE/SE CE/SE
reference 0.6002 reference

0.6

0.5998
h+b

1
φ

0.5996

0.5994

0.5992

0.9995 0.599
−0.5 0 0.5 1 −0.5 0 0.5 1
x−axis x−axis
h+b at t=2 φ at t=2
1.0005 0.6004
CE/SE CE/SE
reference 0.6002 reference

0.6

0.5998
h+b

1
φ

0.5996

0.5994

0.5992

0.9995 0.599
−0.5 0 0.5 1 −0.5 0 0.5 1
x−axis x−axis

67
 h+b at t=4 φ at t=4
1.0005 0.6004
CE/SE CE/SE
reference 0.6002 reference

0.6

0.5998
h+b

φ
0.5996

0.5994

0.5992

0.9995 0.599
−0.5 0 0.5 1 −0.5 0 0.5 1
x−axis x−axis
h+b at t=6.25 φ at t=6.25
1.0005 0.6004
CE/SE CE/SE
reference 0.6002 reference

0.6

0.5998
h+b

1
φ

0.5996

0.5994

0.5992

0.9995 0.599
−0.5 0 0.5 1 −0.5 0 0.5 1
x−axis x−axis

Figure 4.14: Problem 11 : Perturbation of steady-state at rest (h̃ = φ̃ = 10−3 ).

68
 2.1
0.43
2.05

volume fraction
0.42
2
h

0.41
1.95 0.4
1.9 0.39
2 2
2 2
1 1
1 1
y−axis 0 0 y−axis 0 0
x−axis x−axis

1.005
1.005

1
us

1
s
v

0.995 0.995
2 2
2 2
1 1
1 1
x−velocity 0 0 x−axis y−velocity 0 0 x−axis

1.6
1.5
1.4

1.2
uf

1
vf

0.8 0.5
2 2
2 2
1 1
1 1
x−velocity 0 0 x−axis y−velocity 0 0 x−axis

Figure 4.15: Problem 10: Results on 100 × 100 grid points at time t = 0.5.

69
Chapter 5

Conclusion

The space-time CESE scheme was implemented to solve one and two-dimensional SPSF
and TPSF shallow flow models. For authentication, the results of the suggested numerical
scheme are compared with those from the high resolution central scheme. Several case
studies were carried out to check the performance of our scheme and to observe the well
balancing property of our proposed scheme. The results of the suggested numerical scheme
were compared with those from the high resolution central scheme. A close agreement in the
numerical results of both the schemes was observed. However we found that CESE scheme
gives better resolution of the sharp discontinuities as compared to the central schemes and
have comparable accuracy with other schemes available in the literature. We also observed
that our proposed scheme is well balanced and it achieves the steady state conditions.

70
Bibliography

[1] J.J. Dronkers, Tidal computations in rivers and coastal waters, North-Holland Pub.
Co., Amsterdam 1964.

[2] G.A. Corby, Laplace’s tidal equations, an application of solutions for negative depth,
Wiley, 2006.

[3] Charney et al., On the scale of atmospheric motions, (1948)

[4] H.M. Cekirge, R.W. Lardner and R.J. Fraga, Adaptation of the solution of the two-
dimentional tidal equations using the method of characteristics to wind induced cur-
rents and storm surges. Computers and Mathematics with applications, 12(A) (1986)
1081-1090.

[5] S.K. Dube, P.C. Sinha, and G.D. Roy. The numerical simulation of storm surges along
the Bangla Desh coast. Transanctions of the American mathematical society, 9 (1985)
121-133.

[6] A. Chertock, A. Kurganov and G. Petrova. Finite volume particle methods for models
of transport of pollutant in shallow water. Journal of Scientific Computing, 27 (2006)
189-199.

[7] L. Fraccarollo and H. Capart. Riemann wave description of erosional dam-break flows.
Journal of Fluid Mecahnics, 461 (2002) 183-228.

[8] L. Fraccarollo, H. Capart and Y. Zech. A Gudonov method for the computation of
erosional shallow water transients. International Journal for Numerical Methods in
Fluids, 41 (2003) 951-976

71
 [9] C. Synolakis, E. Okal and E. bernard. The megatsunami of December 26,2004. Na-
tional Academy of Engineering Publications,2005.

[10] A.B. Kennedy, Q. Chen, J.T. Kerby and R.A. Dalrymple. Boussinesq modeling of
wave transformation, breaking and runup. Journal of Waterway Port Coastal and
Ocean Engineering, 126(1) (2000) 39-47.

[11] P. Watts, S.T. Grilli, J.T. Kirby, G.J. Fryer and D.R. Tappin. Landslide tsunami
case studies using a Boussinesq model and a fully nonlinear tsunami generation
model.Natural Hazars and Earth System Sciences, 3(5) (2003) 391-402.

[12] C. Zoppou and S. Roberts. Catastrophic collapse of water supply reservoirs in urban
areas. Journal of Hydraulic Engineering, 125(7) (July 1999) 686-695.

[13] H.J.M. Ogink, J.G. Grijsen and A.J.H. Wijbenga. Aspects of flood level computations.
In International Symposium on Flood Frequency and Risk Analysis, USA, 1986.

[14] N. Goutal and F. Maurel. A finite volume solver for 1-D shallow water equations
applied to an actual river. International Journal for Numerical Methods in Fluids, 38
(2002) 1-19.

[15] J.B. Nixon and B.J. Noye. Prawn larvae dispersion modeling. In R.L.May and
A.K.Easton, editors, Computational Techniques and Applications. Computational
Mathematics Group, Australian Mathematical Society, World Scientific, 1996.

[16] K. Hutter and S.B. Savage, Motion of a finite mass of granular material down a rough
incline. J. Fluid. Mech. 199 (1989) 177-215.

[17] F. Bouchut, M. Pelanti, A. Mangeney, A Roe-type scheme for two-phase shallow

granular flows over variable topography, ESAIM-Math. Model. Numer. 42 (2008) 851-
885.

[18] R.M. Iverson, The physics of debris flows, Rev. Geophys. 35 (1997) 245-296.

[19] A. Harten, P.D. Lax, B. van Leer, On upstream differencing and Godunov-type
schemes for hyperbolic conservation laws, SIAM Review 25 (1983) 35-61.

72
[20] E.F. Toro, Shock-capturing methods for free-surface shallow flows, Wiley, Chichester
2001.

[21] C. Berthon, F. Marche, A positive preserving high order V.F.Roe scheme: A class of
Relaxation schemes, SIAM J.Sci. Comp. 30 (2008) 2587-2612

[22] M. Pelanti, A. Mangeney, F. Bouchut, A Riemann solver for SPSF and TPSF models
based on relexation, Comp. Phys. 230 (2011) 515-550

[23] S.C. Chang, The method of space time CESE- New approach for solving Navier Stokes
and Euler equations, J. Comput. Phys. 119 (1995) 295-234.

[24] C.Y. Chow, X.Y. Wang, S.C. Chang, New developments in the method of space-
time CESE-Applications to 2-dimensional time-marching problems, NASA TM (1994)
106758.

[25] C.Y. Chow, X.Y. Wang, S.C. Chang, The space-time CESE. A new high resolution
and genuinely multidimensional paradigm for solving conservation laws, J. Comput.
Phys. 156 (1999) 89-136.

[26] M. Liu. J.B. Wang, K.-Q Wu, The direct aero-acoustics simulation of flow around a
square cylinder using the CE/SE scheme, J. Algorithms Comput. Technol. 1 (2009)
525-537.

[27] S.C. Chang, X.Y. Wang, W.M. To, Application of the space-time CESE method to
1-D convection-diffusion Problems, J. Comput. Phys. 165 (2000) 189-215.

[28] K.B.M.Q. Zaman, C.Y. Loh, Numerical investigation of transonic resonance with
convergent-divergent nozzle, AIAA J. 40 (2002) 2393-2401.

[29] C.Y. Loh, L.S. Hultgren, S.C. Chang, Wave computation in incompresible flow using
the space-time conservatio element and solution element method, AIAA J. 39 (2001)
794-801.

[30] C.Y. Loh, L.S. Hultgren, S.C. Chang, P.C.E. Jorgenson, Noise computation of a shock-
containing supersonic axisymmetric jet by the CESE Method, AIAA Paper, 2000.

73
[31] S. Mudasser, S. Qamar, On application of a variant CESE method for solving 2-D
ideal MHD equations, Appl. Num. Math. 60 (2010) 587-606.

[32] A. Mangeney, M. Pelanti, J.P. Vilotte, F. Bouchut, Numerical modeling of two-

phase gravitational granular flows with bottom topography, in: S. Benzoi-Gavage, D.
Serre (Eds.), Hyperbolic Problems: Theory, Numerics and Applications, Proceedings
of Eleventh International Conference on Hyperbolic Problems (Lyon, France, 2006),
Springer, 2008, 825-832.

[33] A. Kurganov, G. Petrova, Second order well-balanced positivity preserving central-

upwind scheme for the Saint-Venant system. Commun. Math. Sci. 5 (2007) 133-160.

[34] Kyle T. Mandli, Finite Volume methods for the multilayer shallow water equations
with applications to storm surges, 2011.

[35] John Jakeman, On numerical solutions of the shallow water wave equations, October
2006.

[36] Z.C. Zhang, S.T. Yu, S.C. Chang, A Space-time conservation element and solution
element method for solving the two- and three-dimensional unsteady Euler equations
using quadrilateral and hexahedeal meshes. J. Comput. Phys., 175 (2002) 168-199.

[37] G. Warnecke, S. Qamar, Application of space-time CESE method to shallow water

magnetohydrodynamic squations. J. Comput. Appl. Math. 196 (2005) 132-149.

[38] C.B. Vreugdenhil, 2-layer shallow -water flow in 2 dimensions, a numerical study. J.
Comput. Phys. 33 (1979) 169Ű184.

[39] R.J. LeVeque, Balancing source terms and flux gradients in high-resolution Godunov
methods: The quasi-steady wavepropagation algorithm. J. Comput. Phys. 146 (1998)
346Ű365.

[40] R. Abgrall, S. Karni, 2-layer shallow water system: a relaxation approach. SIAM J.
Sci. Comput. 31 (2009) 1603-1627.

74
[41] K. Xu. A well-balanced gas-kinetic scheme for the shallow-water equations with source
terms. J. Comput. Phys. 178 (2002) 533-562.

[42] M.J. Castro, J. Mact’ýas and C. Part’es, A Q-scheme for a class of systems of coupled
conservation laws with source term. Application to a two-layer 1-D shallow water
system. ESAIM: M2AN 35 (2001) 107Ű127.

[43] M.J. Castro, J.A. Garct’ýa Rodrt’ýguez, J.M. Gonzt’alez-Vida, J. Mact’ýas, C.

Part’es,and M.E. Vt’azquez-Cendt’on, Numerical simulation of 2 layer shallow water
flows through channels with irregular geometry. J. Comput. Phys. 195 (2004) 202Ű235.

[44] F. Bouchut, Non-linear stability of finite volume methods for hyperbolic conservation
laws and well-balanced schemes for sources, Birkhauser-verlag, 2004.

[45] S. Jin, A steady-state capturing method for hyperbolic system with geometrical source
terms, Math. Model. Anal. 35 (2001) 631-646.

[46] S. Qamar, S. Mudasser. A Kinetic flux-vector splitting method for the shallow water
magnetohydrodynamics. J.Comp. Phys. 181 (2010) 1109-1122

75