NUMERICAL MODELLING

OF
TSUNAMI WAVE EQUATIONS

A
thesis

submitted in partial fulfillment of the

requirements for the award of the degree
of

MASTER OF SCIENCE
IN
MATHEMATICS

submitted by
HARI SHANKAR SHAW
Roll No-412MA2085

under the supervision

of
Prof. SNEHASHISH CHAKRAVERTY

DEPARTMENT OF MATHEMATICS
NIT ROURKELA
ROURKELA– 769 008

MAY 2014
 NATIONAL INSTITUTE OF TECHNOLOGY
ROURKELA

DECLARATION

I hereby declare that the work which is being presented in the report entitled “ NU-
MERICAL MODELLING OF TSUNAMI WAVE EQUATIONS ” in par-
tial fulfillment of the requirement for the award of the degree of Master of Science,
submitted in the Department of Mathematics, National Institute of Technology,
Rourkela is an authentic record of my own work carried out under the supervision
of Prof. S. Chakraverty.

The matter embodied in this has not been submitted by me for the award of any
other degree.

May, 2014

(Hari Shankar Shaw)

 CERTIFICATE

This is to certify that the project report entitled “ NUMERICAL MODELLING

OF TSUNAMI WAVE EQUATIONS ” submitted by Hari Shankar Shaw
to the National Institute of Technology Rourkela, Odisha for the partial fulfillment
of requirements for the degree of master of science in Mathematics is a bonafide
record of review work carried out by him under my supervision and guidance. The
contents of this report, in full or in parts, have not been submitted to any other
institute or university for the award of any degree or diploma.

May, 2014
(Prof. S. CHAKRAVERTY)
Head of Department
Department of Mathematics
NIT Rourkela
(Hari Shankar Shaw)

ii
ACKNOWLEDGEMENTS

I am extremely indebted to be involved in a challenging research project like “Nu-

merical Modelling Of Tsunami Wave Equations”. This project increased my thinking
and understanding capability as I started the project from scratch.
I would like to express my greatest gratitude and respect to my supervisor Prof.
Snehashish Chakraverty, for his valuable suggestions, excellent guidance and endless
support. He is not only a wonderful supervisor but also a genuine person. I consider
myself extremely lucky to be able to work under guidance. Actually he is one of such
genuine person for whom my words are not enough to express.
I am very much grateful to Prof. Sunil Kumar Sarangi, Director, National In-
stitute of Technology, Rourkela for providing excellent facilities in the institute for
carrying out research.
I would like to express my sincere thanks to PhD scholars especially Sukant
Nayak and Karan Kr. Pradhan for their help and precious suggestions to perform
the project work. I am very much thankful to them for giving his valuable time for
me.
I would like to express my thanks to all the faculty members, all my classmates,
all staffs of Mathematics department for making my stay in N.I.T. Rourkela a pleas-
ant and memorable experience and also giving me absolute working environment
where I unleashed my potential.
I owe a gratitude to God and my parents for their blessings and inspiration.

iii
ABSTRACT
This report investigates the modelling of tsunami wave using one dimensional shal-
low water equations (SWEs) by numerical methods namely finite difference method
(FDM) and finite volume method (FVM). We have used one dimensional SWEs
to model the water wave propagation i.e. we study the variation of water surface
elevation with finite distance. We obtained the SWEs from Euler’s equation of mass
and momentum assuming a long wave approximation. First of all we approximate
the SWEs using FDM and then by FVM for showing the behaviour of water surface
elevation with distance. After approximating the SWEs using both the numerical
method, results have been shown using different schemes viz. FDM as well as FVM.
Moreover, in actual practice, we may have incomplete information about the vari-
ables being a result of errors in modelling, observations, or by applying different
initial as well as boundary conditions etc. Rather than the particular value of water
surface we may have only the bounds of the values. These bounds may be given in
term of interval. Thus we have developed interval finite volume method (IFVM) also
for approximating one dimensional SWEs to model tsunami wave with uncertain (in-
terval) parameter. Next, numerical results have been shown using IFVM. Then a
comparision study has been investigated to compare the results of both the method
i.e FDM and FVM. Finally all computed results are shown in terms of tables and
plots.

iv
Contents
1 Intoduction 1

2 One dimensional shallow water equations (SWEs) 4

3 Finite difference method (FDM) for solving one dimensional SWEs 5

3.1 Intoduction to finite difference method (FDM) . . . . . . . . . . . . . 5
3.2 Different schemes of FDM . . . . . . . . . . . . . . . . . . . . . . . . 6
3.3 Explicit scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
3.4 Semi-implicit scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
3.5 Implicit scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.6 Tsunami wave approximation . . . . . . . . . . . . . . . . . . . . . . 8

4 Finite volume method (FVM) for solving one dimensional SWEs 9

4.1 Introduction to finite volume method (FVM) . . . . . . . . . . . . . . 9
4.2 Different schemes of FVM . . . . . . . . . . . . . . . . . . . . . . . . 11
4.3 Upwind interpolation . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
4.4 Central difference (CD) interpolation . . . . . . . . . . . . . . . . . . 13
4.5 Tsunami wave approximation . . . . . . . . . . . . . . . . . . . . . . 14
4.6 Interval finite volume method (IFVM) . . . . . . . . . . . . . . . . . 16

5 Numerical results 17
5.1 Numerical results using FDM . . . . . . . . . . . . . . . . . . . . . . 17
5.2 Numerical results using FVM and IFVM . . . . . . . . . . . . . . . . 20

6 Comparison between FDM and FVM, conclusion and future direc-

tions 23

References 26

v
1 Intoduction
Tsunamis are generated by the movement of sea bottom due to long waves of earth-
quakes. The impulsive sea floor movement in the earthquakes region causes the
water surface region instantaneously as discussed in [1]. The sudden gain in po-
tential energy converts to kinetic energy by the gravitational force which serves as
the restoring force of the system. Generally, tsunamis are treated as shallow water
waves.
Tsunami is a Japanese word that is the combination of two words :“ tsu ” means
harbor and “ nami ” means wave. Therefore, tsunami literally means “harbor
waves”. The word was originally created to describe large amplitude oscillations
in a harbor under the resonance condition given in [2]. The most common cause of
tsunami are under sea earthquakes.
In this report we have used shallow water equations ( SWEs ) to model water
wave propagation in one dimension. We obtain the SWEs from Euler’s equation of
mass and momentum considering a long wave approximations given in [1], [3]. SWEs
state the propagation of water waves whose wave length is much longer than the
depth of water. Therefore, we have modeled tsunami wave using SWEs.

Shallow Water Equations (SWEs)

Shallow Water Equations ( SWEs ) are a system of hyperbolic partial differen-

tial equations ( PDEs ) governing the flow of fluid in the rivers, channels, oceans
and costal regions. We have investigated SWEs from mass and energy conservation
principle expressed in the Navier-Stokes equations. SWEs give the idea about the
flow of water waves, especially those water wave whose wave length is much longer
than the depth ( basin ) of water.
The wavelength of tsunami waves are far longer than the normal waves. A tsunami
wave initially resembles a rapidly rising tides for this reason they are often referred
to as tidal waves. The average depth [3] of ocean nearly 5 Km, which is compared
with the wavelength long waves or tsunamis, which may exceed 100 Km.
Although the impact of tsunamis is limited to coastal areas, their destructive
power can be enormous and they can affect entire ocean basins; in 2004 Indian
Ocean tsunami was among the deadliest natural disasters in human history with
over 230,000 people killed in 14 countries bordering the Indian ocean. One can get
SWEs by neglecting the bottom friction and assuming long wave approximations
from the Euler equations of mass and momentum.
SWEs are used to model Dam Beaks, strom surges, solute transport, river flows, eco-
nomical model etc given in [3]. The problem with SWEs is, they are difficult to model
in dry areas where water is not present. SWEs are only defined in wet regions. Thus
for these type of equations we need actually to deal with moving boundary problems.

How do SWEs arise ?

SWEs are investigated by Navier-Stokes (N-S) equations, which describe the mo-
tion of fluid. Also, the N-S equations are derived from the equations of conservation

1
of mass and linear momentum discussed by Imamura [1]. So, from the translation
motion of fluid element and neglecting the vertical acceleration, the equations of
mass conservation and momentum in one dimensional problem are described as fol-
lows [1]  ∂u ∂w
 ∂x + ∂z = 0






∂u
∂t
+ u ∂u
∂x
+ w ∂u
∂z
+ (1/ρ) ∂P
∂x
=0 (1)






 ∂w
∂t
+ u ∂w
∂x
+ w ∂w
∂z
+ g + (1/ρ) ∂P
∂z
=0

where, x is the horizontal axis and z is the vertical axis, t is time, η is the water
surface elevation, h is basin depth, u, w are the velocities of fluid in the x and z
directions respectively, and g is the acceleration due to gravity. Finally after taking
the dynamics and kinetic conditions at surface and bottom [3], we get the Shallow
Water Equations (SWEs).
Mathematical modeling plays a vital role in the area of tsunami science, such as
in the area of scientific studies for tsunami propagation and initiation. This thesis
investigates the numerical solution of one dimensional SWEs using numerical meth-
ods namely Finite Difference Method (FDM) and Finite Volume Method (FVM).

FVM is one of the most useful method for modeling the SWEs, long waves, ra-
diative transfer, etc. FVM is widely used in engineering, fluid mechanics, petroleum
engineering, computational fluid dynamics, heat and mass transfer, etc. The most
important feature of this method is that numerical flux is conserved from one dis-
cretization cell to its neighbor. This feature makes FVM quite attractive for modeling
problems in fluid mechanics, heat transfer and semi conductor device simulation.
FVM is a method of representing and evaluating the partial differential equa-
tions to algebraic equations. In this method we calculate the values at discrete
places on meshed geometry as in finite difference method (FDM) or finite element
method (FEM). Finite (control) refers to small volume surrounding each node point
on a mesh. Then we used interval finite volume method (IFVM) to the SWEs and
have investigated the variation of Tsunami wave.

The shallow water equations introduced in [6] is very commonly used for nu-
merical solution. Modelling of tsunami wave has been solved using linear Leap-frog
method by Imamura, Yalciner [1]. A detailed study of SWEs is given in IUGG/IOC
time project [3] for numerical methods in tsunami simulation. A numerical modelling
on one dimensional and two dimensional SWEs using FDM discussed in Junbo Park
[4] and added the numerical simulation of wave propagation. A lot of research work
on one and two dimensional convection-diffusion problems, [8] has been solved using
finite volume scheme by Versteeg and Malalasekera. Cebeci, et al. [9] proposed real
life problems using FVM.

In this report our main aim is to develop an efficient numerical method for solving
SWEs to model Tsunami wave. Generally, the values of variables or properties are

2
taken as crisp but in actual case the accurate (crisp) values may not be obtained.
To overcome the vagueness we use interval in place of crisp values. So, next aim is
to study the interval finite volume method (IFVM). The IFVM has been developed
here to study the variation of tsunami waves.
In this report we first discuss the introduction of tsunami wave and shallow water
waves and their origins in section 1 and 2. A detailed study of one dimensional SWEs
using different schemes of FDM has been done in section 3. Then we have solved one
dimensional SWEs using FVM with different schemes of interpolation and variation
for FVM and IFVM in section 4. Section 5 deals with the numerical results for one
dimensional SWEs using FDM and FVM . A comparative numerical results using
FDM and FVM has been investigated in section 6. Also conclusion and future work
has been included in section 6 and finally references are cited.

3
2 One dimensional shallow water equations (SWEs)
The shallow water Eqs. in one dimension is given by [1]
∂η ∂M


 + =0 (mass conservation law)
 ∂t
 ∂x
(2)
 ∂M + gD ∂η = 0 (momentum conservation law)



∂t ∂t
where,

η = Water surface elevation

M = Discharge flux in the positive x-direction

g = Acceleration due to gravity

D = Total thickness of water

h = Basin depth of water

Thus, D = η + h

Fig. 1 shows the behaviour of one dimensional shallow water equations, where
vertical axis represents water surface elevation and horizontal axis represents dis-
tance.

Figure 1: 1-D shallow water model

The Eqs. in (2) are coupled first order partial differential equations which can
be uncoupled to produce two second order partial differential equations [4] which
are given as follows
 2
∂ M ∂ 2M 1 ∂M ∂M
= g(η + h) −


2 2

 ∂t ∂x (η + h) ∂x ∂t


(3)
 2 2
 2

 ∂ η ∂ η ∂η ∂h ∂η
 2 = g(η + h) 2 + g +g


∂t ∂x ∂x ∂x ∂x

4
3 Finite difference method (FDM) for solving one
dimensional SWEs
In this report we have used FDM for solution of one dimensional SWEs.

3.1 Intoduction to finite difference method (FDM)

Finite-difference methods are numerical methods for approximating the solutions to
differential equations using finite difference equations to approximate derivatives.

From Taylors Series expansion, we have

∂Φ 1 ∂ 2Φ 2 1 ∂ nΦ
Φ(x + ∆x) = Φ(x) + (∆x) + (∆x) +, ..., + (∆x)n . (4)
∂x 2! ∂x2 n! ∂xn

The grid generation is well known in FDM and so has been depicted from Fig. 2

Figure 2: Grid generation in FDM

In view of Fig. 2 and Eq. (4) we can write

 j  2 j+1  3 j
j+1 j ∂Φ ∂ Φ 2 ∂ Φ
Φi = Φi + ∆x + 2
(1/2!)(∆x) + (1/3!)(∆x)3 + ...
∂x i ∂x i ∂x3 i

j j+1 j
(Φj+1 − Φji ) ∂ 2Φ ∂ 3Φ
  
∂Φ
⇒ = i − (1/2!)∆x − (1/3!)(∆x)2 . (5)
∂x i ∆x ∂x2 i ∂x3 i

As such forward difference approximation may be written as

Φji+1 − Φji
 j

∂Φ
= + (∆x) (6)
∂x i ∆x
Similarly we have, backward difference approximation as

Φji − Φji−1
 j

∂Φ
= + (∆x) (7)
∂x i ∆x

5
 Finally the central difference approximation is
j j+1 j
∂ 2Φ ∂ 3Φ
  
∂Φ
Φj+1
i = Φji + ∆x + (1/2!)(∆x) + 2
(1/3!)(∆x)3 + ... (8)
∂x i ∂x2 i ∂x3 i

and
j j+1 j
∂ 2Φ ∂ 3Φ
  
∂Φ
Φji−1 = Φji − ∆x + (1/2!)(∆x) − 2
(1/3!)(∆x)3 + ... (9)
∂x i ∂x2 i ∂x3 i

Subtracting (9) from (8), we get

 j
∂Φ
Φji+1 − Φji−1 =2 ∆x + (∆x)2
∂x i

Φji+1 − Φji−1
 j

∂Φ
⇒ = + (∆x)2 (10)
∂x i 2∆x

3.2 Different schemes of FDM

We have used different schemes of FDM such as Explicit, Semi-implicit, and Implicit
schemes for numerical solutions of one dimensional shallow water equations.

3.3 Explicit scheme

To solve the SWEs in one dimension we discretize the first Eq. (2) with respect to
both space and time. We approximate the time derivative by forward difference and
space derivative by central difference, thus we have
∂η ∂M
+ = 0
∂t ∂x
" # " #
j j
ηij+1 − ηij Mi+1 − Mi−1
⇒ + =0
∆t 2 ∆x
∆t  j
⇒ ηij+1 − ηij = − j
  
Mi+1 − Mi−1
2 ∆x

⇒ ηij+1 = ηij − (c/2) Mi+1

 j j

− Mi−1 (11)

Also from second Eq. of (2), we have

∂M ∂η
+ gD = 0
∂t ∂x

6
 " # " #
j j
Mij+1 − Mij ηi+1 − ηi−1
⇒ + gD =0
∆t 2 ∆x

⇒ Mij+1 = Mij − (1/2)c g (ηij + h) ηi+1

 j j

− ηi−1 (12)

Initially we are taking the basin depth to be zero. i.e. h = 0. So from the above the
equation we get

⇒ Mij+1 = Mij − (1/2)c g (ηij ) ηi+1

 j j

− ηi−1 (13)

∆t
where, c = , which is the ratio between time step and spatial step.
∆x

3.4 Semi-implicit scheme

Using this scheme of FDM we can approximate [4] the first Eq. given in (2) as follows

∂η ∂M
+ = 0
∂t ∂x
" # " #
j
ηij+1 − ηi−1 j
Mi+1 j
− Mi−1
⇒ + =0
∆t 2 ∆x
∆t  j
⇒ ηij+1 − ηij = − j
  
Mi+1 − Mi−1
2 ∆x

⇒ ηij+1 = ηij − (c/2) Mi+1

j j


− Mi−1 (14)
Now applying Crank-Nicolson approximation to the Eq. (14) we have
j j

 j j j+1 j+1

Mi+1 − Mi−1 = 1/2 (Mi+1 − Mi−1 ) + (Mi+1 − Mi−1 )

⇒ ηij+1 = ηij − (c/4)[Mi+1

j j
− Mi−1 j+1
] − (c/4)[Mi+1 j+1
− Mi−1 ] (15)
Again from the second Eq. of (2), we have
∂M ∂η
+g D = 0
∂t ∂x
" # " #
j j
Mij+1 − Mij ηi+1 − ηi−1
⇒ +g D =0
∆t 2 ∆x
⇒ Mij+1 = Mij − (1/2)c g ηij + h ηi+1
  j j

− ηi−1

⇒ Mij+1 = Mij − (1/2)c g ηij ηi+1

  j j

− ηi−1 (16)

7
3.5 Implicit scheme
We now implement another scheme of FDM known as implicit method [4]on the
one dimensional SWEs. In this method we have used Crank-Nicolson approxima-
tion, which is the average of the central differences about the point (i, j) and (i, j +
1) . This method has a stable solution for any value of ∆t and ∆x . So, after applying
this method in the first Eq. of (2) and solving we get
∂η ∂M
+ = 0
∂t ∂x
j j
ηij+1 − ηij Mi+1 − Mi−1
⇒ + =0
∆t 2 ∆x
∆t
⇒ ηij+1 − ηij = − j j



Mi+1 − Mi−1
2 ∆x
⇒ ηij+1 = ηij − (c/2) Mi+1j j


− Mi−1 (17)
Again applying Crank-Nicolson approximation to the above Eq. we have
j j

 j j j+1 j+1

Mi+1 − Mi−1 = 1/2 (Mi+1 − Mi−1 ) + (Mi+1 − Mi−1 )
⇒ ηij+1 = ηij − (c/4)(Mi+1
j j
− Mi−1 j+1
) − (c/4)(Mi+1 j+1
− Mi−1 )

⇒ ηij+1 + (c/4)(Mi+1
j+1 j+1
− Mi−1 ) = ηij − (c/4)(Mi+1
j j
− Mi−1 ) (18)
From the second Eq. of (2) now we have
∂M ∂η
+g D = 0
∂t ∂x
" # " #
j j
Mij+1 − Mij ηi+1 − ηi−1
⇒ + gD =0
∆t 2 ∆x
⇒ Mij+1 = Mij − (1/2)cg ηij + h ηi+1

j j


− ηi−1

⇒ Mij+1 +(1/4)cgh ηi+1

 j+1 j+1
= Mij −(1/2)cg ηij ηi+1

 j j
  j j

− ηi−1 − ηi−1 −(1/4)cgh ηi+1 − ηi−1
(19)
Left hand sides of Eqs. (18) and (19) are the terms of η and M at time
step j + 1. The right hand sides are at time step j. So, we can solve the above
system of Eqs. by any well known method.

3.6 Tsunami wave approximation

In previous sections, we have assumed the basin depth of water to be zero. i.e. h =
0 . Now we have approximated the basin depth h by a hyperbolic tangent [4] as
given below

(x − 70)
h(x) = 50 − 45 tanh [ ] where, 0 m 6 x 6 100 m(20)
8

8
So, from Eq. (20) we can find the maximum and minimum value of h as 95 m and
5 m respectively.
Thus, for the explicit scheme case i.e. from Eq. (13) we have

⇒ Mij+1 = Mij − (1/2)c g (ηij ) ηi+1

 j j

− ηi−1
⇒ Mij+1 = Mij − (1/2)c g (ηij + h) ηi+1
 j j

− ηi−1 (21)

and for the implicit scheme case i.e. from Eq. (19) we have

⇒ Mij+1 = Mij − (1/2)cg ηij + h ηi+1

j j


− ηi−1
⇒ Mij+1 = Mij − (1/2)c g ηij ηi+1

j j

j j


− ηi−1 − 1/2c g h ηi+1 − ηi−1

⇒ Mij+1 = Mij −(1/2)c g ηij j j

 j j j+1 j+1


ηi+1 − ηi−1 −1/4c g h (ηi+1 − ηi−1 ) + (ηi+1 − ηi−1 ]

⇒ Mij+1 +1/4c g h ηi+1

j+1 j+1
= Mij −(1/2)c g ηij ηi+1


j j

j j


− ηi−1 − ηi−1 −1/4c g h ηi+1 − ηi−1 (22)

4 Finite volume method (FVM) for solving one

dimensional SWEs
In this report we have used another numerical method known as finite volume
method (FVM) for solution of one dimensional SWEs.

4.1 Introduction to finite volume method (FVM)

FVM is a method of representing the partial differential equations into algebraic
equations. In this method, we calculate the values at discrete places or points as
in finite difference or finite element methods. Finite (control) refers here as small
volume surrounding each node point on a mesh.
In this method, we integrate the given partial differential equations that con-
tains a divergence term which can be converted to surface integral using divergence
theorem. FVM is based on discretization of the integral forms of the conservation
equations. Discretization is applied directly to the integral equations for small con-
trol (finite) volumes as shown in the Fig. 3 [7]

9
 Figure 3: One dimensional control (finite) volume in FVM

In this method, instead of discretizing first, we start with the integral form of
the equations.

Below we write the steps involved in FVM for solving one dimensional SWEs.

Step 1-Grid generation in FVM:

The first step in the finite volume method is grid generation i.e. by dividing the
domain into discrete control volumes. The boundaries of control volumes are posi-
tioned mid-way between adjacent nodes. Thus each node is surrounded by a control
volume or cell. Generally, it is better to set up control volumes near the edge of the
domain in such a way that the physical boundaries coincide with the control volume
boundaries. Consider a control volume whose nodal point is P and the neighbouring
the nodes to the west and east, are defined as W and E respectively. The west face
of the control volume is referred by w and east face by e. The distance from W to
P , and P to E are given by δxW P and δxP E respectively as shown in the Fig. 4. Also
the distance from w to P , P to e are given by δxwP and δxP e respectively. The
width of control volume [8] from w to e is denoted as ∆x = δxwe as shown in the
fig (4).

Figure 4: One dimensional grid for FVM [8]

Step 2-Discretization Concept:

10
 The key step of FVM is the integration of the governing equation over a fi-
nite ( control ) volume to yield a discretised equation at its nodal point P .

Step 3- Solution of the problem:

Discretised equation must be set at each of the nodal points in order to solve
a problem.The resulting system of linear algebraic equations is then solved by any
well known numerical method.

4.2 Different schemes of FVM

We have used two schemes of FVM namely upwind interpolation (UI) and central
difference (CD) interpolation method for solving one dimensional shallow water
equations .

4.3 Upwind interpolation

It is the simplest way for approximating the partial differential equations by FVM. Here, we
use the value at a neighboring grid point. We have taken the nearest upwind ( up-
stream ) grid point. Thus, taking the integral form of first Eq. of (2) into small
control ( finite ) volumes, and then discretized in the nearest grid point. The volume
integral is converted into surface integral using divergence theorem. Thus, from first
Eq. (2) we have
∂η ∂M
+ = 0 (23)
∂t ∂x

Integrating over control (finite) volume, we get

d xi+1/2
Z
ηdx + Me − Mw (24)
dt xi−1/2

Now we approximate the above Eq. [7]

d xi+1/2
Z
ηdx ≈ ηP ∆x (25)
dt xi−1/2

At the cell boundary i.e. in the east face e = xi+1/2 , the normal n is in the positive
direction, so

Me ≈ MP

Again at the cell boundary i.e. at the west face w = xi−1/2 the normal is in the
negative direction. So, taking the value at the nearest grid point in the west of the
cell we have

Mw ≈ MW

11
 So the finite volume approximation of the above Eq. (24) is
d
(ηP ∆x) + MP − MW = 0
dt
" # " #
j+1 j j j
η − ηi Mi − Mi−1
⇒ i + =0
∆t ∆x

⇒ ∆x ηij+1 − ηij = ∆t Mij − Mi−1 j

⇒ ηij+1 = ηij − c Mij − Mi−1

j
 
(26)
∆t
where, c = , which is the ratio of time step and spatial step.
∆x
Again from the second Eq. of (2) we have
∂M ∂η
+ gD =0 (27)
∂t ∂x
Integrating over control (finite) volume, we get
d xi+1/2
Z Z xi+1/2
∂η
M dx + g(η + h) dx = 0 (28)
dt xi−1/2 xi−1/2 ∂t

where, g is the acceleration due to gravity, and h is the basin depth of water. Taking
h to be zero, the Eq. becomes
d xi+1/2
Z Z xi+1/2
∂η
M dx + g(η) dx = 0
dt xi−1/2 xi−1/2 ∂t

d xi+1/2
Z Z xi+1/2
∂η
⇒ M dx + g η dx = 0
dt xi−1/2 xi−1/2 ∂x
d g
⇒ (MP ∆x) + [(ηe + ηw ) (ηe − ηw )] = 0
dt 2
" #
Mij+1 − Mij
 
g (ηp + ηw )(ηp − ηw )
⇒ + =0
∆t 2 ∆x
" # " #
j j j j
Mij+1 − Mij g (ηi + ηi−1 )(ηi − ηi−1 )
⇒ + =0
∆t 2 ∆x
∆t g  j
⇒ Mij+1 − Mij = j
)(ηij − ηi−1
j
  
(ηi + ηi−1 ) =0
∆x 2
 cg j
⇒ Mij+1 − Mij = j
)(ηij − ηi−1
j
 
(ηi + ηi−1 ) =0
2
cg  j
⇒ Mij+1 = Mij − j
)(ηij − ηi−1
j

(ηi + ηi−1 ) (29)
2
∆t
where, c =
∆x

12
4.4 Central difference (CD) interpolation
This approximation is based on central difference interpolation between two neigh-
boring grid points. Using this interpolation we solve the first Eq. of (2) in the
following way
Integrating over control (finite) volume of the grid points as shown in Fig. 4 we have
d xi+1/2
Z Z xi+1/2
∂M
ηdx + dx = 0 (30)
dt xi−1/2 xi−1/2 ∂x

Approximating we have [7]

Z xi+1/2
d
ηdx ≈ ηP ∆x (31)
dt xi−1/2

Now, the x derivatives are approximated by central difference interpolation. Accord-

ingly we have the following
" #
ηij+1 − ηij
∆x + (M )e − (M )w = 0
∆t
" #
ηij+1 − ηij
⇒ ∆x + [(M )E − (M )W ] = 0
∆t
" # 
ηij+1 − ηij

(M )E − (M )W
⇒ + =0
∆t ∆x
" # " #
j j
ηij+1 − ηij Mi+1 − Mi−1
⇒ + =0
∆t ∆x

⇒ ηij+1 = ηij − c Mi+1

 j j

− Mi−1 (32)
Again from the second Eq. of (2) we have
∂M ∂η
+ gD =0 (33)
∂t ∂x

Integrating over control (finite) volume, we get

d xi+1/2
Z Z xi+1/2
∂η
M dx + g(η + h) dx = 0 (34)
dt xi−1/2 xi−1/2 ∂x
Assume h to be zero, the Eq. becomes

Z xi+1/2 Z xi+1/2
d ∂η
M dx + g η ∆x = 0 (35)
dt xi−1/2 xi−1/2 ∂x

13
Also, Z xi+1/2
d
M dx ≈ MP ∆x (36)
dt xi−1/2

The x derivatives are approximated by central difference interpolation. Thus, we

have the following Eqs. as

d
[MP (∆x)] + g I = 0 (37)
dt

where,
Z xi+1/2
∂η
I= η dx (38)
xi−1/2 ∂x

Now, solving I near the neighbouring grid points we have

1
I= [(ηE + ηW )(ηE − ηW )]
2
1 j j j j

⇒I= (ηi+1 + ηi−1 )(ηi+1 − ηi−1 ) (39)
2
Putting the value of I in the Eq. (37) we get

d 1 j j j j

(MP ∆x) + g (ηi+1 + ηi−1 )(ηi+1 − ηi−1 ) =0
dt 2
" #
Mij+1 − Mij g j j j j

⇒ ∆x + (ηi+1 + ηi−1 )(ηi+1 − ηi−1 ) =0
∆t 2
∆t g  j
⇒ Mij+1 − Mij + j j j

(ηi+1 + ηi−1 )(ηi+1 − ηi−1 ) =0
∆x 2
cg  j
⇒ Mi = Mij −
j+1 j

j j


ηi+1 + ηi−1 ηi+1 − ηi−1 (40)
2
The values of M in the left hand side of Eq. (40) represents the values at time
step j + 1 and the right hand side term of Eq. represents the values at time step j.
Thus we can solve Eqs. (32) and (40) for evaluating the values of M and η by an
iterative method.

4.5 Tsunami wave approximation

Initially, in FVM we have assumed the basin depth of water to be zero. i.e. h =
0 . Now we have approximated the basin depth h by a hyperbolic tangent (as taken

14
in FDM)

(x − 70)
h(x) = 50 − 45 tanh [ ] where, 0 m 6 x 6 100 m(41)
8

Thus, from Eq. (28) we have

d xi+1/2
Z Z xi+1/2
∂η
M dx + g (η + h) dx = 0 (42)
dt xi−1/2 xi−1/2 ∂x

Let us assume,
Z xi+1/2
∂η
I= (η + h)
xi−1/2 ∂x
Z xi+1/2 Z xi+1/2
∂η ∂η
⇒I= η + h (43)
xi−1/2 ∂x xi−1/2 ∂x

Now,
Z xi+1/2
∂η 1
η = [(ηe + ηw ) (ηe − ηw )] (44)
xi−1/2 ∂x 2

and,
Z xi+1/2
∂η
h = h (ηe − ηw ) (45)
xi−1/2 ∂x

So,
1
I= [(ηe + ηw ) (ηe − ηw )] + h (ηe − ηw ) (46)
2

Now, putting the value of I in Eq. (42) and solving we get,

d g
(MP ∆x) + [(ηe + ηw ) (ηe − ηw )] + gh (ηe − ηw ) = 0
dt 2
" #
Mij+1 − Mij g
⇒ + [(ηp + ηw )(ηp − ηw ) + gh(ηp − ηw )] = 0
∆t 2∆x
 cg  j
⇒ Mij+1 − Mij + j
)(ηij − ηi−1
j
) + gh ηij − ηi−1
j
   
(ηi + ηi−1 =0
2
cg  j
⇒ Mij+1 = Mij − j
)(ηij − ηi−1
j
) − gh ηij − ηi−1
j
  
(ηi + ηi−1 (47)
2

15
4.6 Interval finite volume method (IFVM)
The uncertain values occurred in practical cases (such as errors in experimental
data, and partial or imperfect knowledge of the parameters) may be handled by
taking the uncertainty as interval sense. So to compute these uncertainties we need
interval arithmetic. Let us consider the uncertain values of a parameter η in interval
form and the same may be written in the following way

[ x , x ] = [ x | x ∈ R, x ≤ x ≤ x ] (48)
where, [ x ] and [ y ] are lower and upper values of the interval respectively. Let us
assume two intervals [ x , x ] and y , y ] then we have
 
[ x ,x ] + y ,y ] = x + y ,x + y ]
 
[ x ,x ] − y ,y ] = x − y ,x − y ] (49)

Upwind Interpolation IFVM:

Now applying the IFVM in Eqs. (26) and (29) we get

[ M , M ]j+1 = [ M , M ]ji − (cg/2) [ η, η ]ji + [ η, η ]ji−1 [ η, η ]ji − [ η, η]ji−1

  
i

[ η, η ]j+1
i = [ η, η ]ji − c[ M , M ]ji − [ M , M ]ji−1 (50)
Rearranging the above equations we get a set of four Eqs. as follows

[ M ]j+1
i = [ M ]ji − (cg/2) [ η ]ji + [ η ]ji−1 [ η ]ji − [ η ]ji−1
[ M ]j+1
i = [ M ]ji − (cg/2) [ η ]ji + [ η ]ji−1 [ η ]ji − [ η ]ji−1
[ η ]j+1
i = [ η ]ji − c [ M ]ji − [ M ]ji−1
[ η ]j+1
i = [ η ]ji − c [ M ]ji − [ M ]ji−1 (51)
Central Difference Interpolation IFVM:

Applying the IFVM the Eqs. (32) and (40) we have

[ M , M ]j+1 = [ M , M ]ji − (cg/2) [ η, η ]ji+1 + [ η, η ]ji−1 [ η, η ]ji+1 − [ η, η]ji−1

  
i

[ η, η ]j+1 = [ η, η ]ji − c M , M ]ji+1 − [ M , M ]ji−1 (52)

 
i

Rearranging the above Eqs. one can get a set of four Eqs. as below

[ M ]j+1 = [ M ]ji − (cg/2) [ η]ji+1 + [ η ]ji−1 [ η ]ji+1 − [ η ]ji−1

  
i

[ M ]j+1 = [ M ]ji − (cg/2) [ η ]ji+1 + [ η ]ji−1 [ η ]ji+1 − [ η ]ji−1

  
i

[ η ]j+1 = [ M ]ji − c [ M ]ji+1 − [ M ]ji−1

 
i

[ η ]j+1 = [ M ]ji − c [ M ]ji+1 − [ M ]ji−1

 
i (53)

16
5 Numerical results
As mentioned earlier we have used finite difference method ( FDM ) and finite volume
method ( FVM ) for solution of one dimensional SWEs. Moreover interval finite
volume method (IFVM) has also been developed to solve the SWEs in uncertain
environment.

5.1 Numerical results using FDM

For one dimensional shallow water equations we have shown the 
numerical
 results
1
for different schemes of FDM. We have used [4] grid size ∆x = m for space
  8
1
and ∆t = s for time step.
3000
The boundary conditions [4] have been taken as

at x = 0 and x = L M = 0.

Also η(0, j) = η(0, j) ; η(L − 1, j) = η(L, j)

The initial conditions are taken as

at t = 0, assuming the initial velocity of water is zero, i.e. water is at rest

position at t = 0 and therefore M = 0. First we take the basin depth of water to
be zero i.e. h = 0. Another assumption is, the maximum and minimum values of
η are 18 m and 20 m respectively.
And the corresponding MATLAB program has been developed to compute the
result for the behavoiur water surface elevation (i.e. η) with distance x is shown in
Fig. 5 at (h = 0, t = 0)

Figure 5: At t = 0s, graph of η with distance x

We have shown the graphs for variation of water surface elevation η with different
value of time t in case of explicit scheme of FDM in Fig.6 to 10 for different value
of t.

17
 Figure 6: At time t = 0.1 s Figure 7: At time t = 0.16667 s

Figure 8: At time t = 0.3333 s Figure 9: At time t = 1 s

Figure 10: At time t = 1.6666 s

Again we have shown the plots for variation of water surface elevation η with
different values of time t in case of implicit scheme of FDM in Fig.11 to 15.

18
 Figure 11: At time t = 0.1 s Figure 12: At time t = .16667 s

Figure 13: At time t = 0.333 s Figure 14: At time t = 1 s

Figure 15: At time t = 1.6667 s

Initially we have assumed the basin depth of water i.e. h = 0. Now from
Eq. (20) we take the minimum and maximum value of basin depth as h = 5 m and h =
95 m respectively. Accordingly we have shown the behaviour of η with distance in
the Fig. 16 and 17 for both the cases.

The behaviour of water surface elevation η with distance x when the basin depths
are 85m 90m and 95m is depicted in Fig. 18.

19
Figure 16: At t = 0.1s, graph of Figure 17: At t = 0.1s, graph of
η with distance x when h = 5m η with distance x when h = 95m

Figure 18: graph for h = 85m, 90m, 95m

5.2 Numerical results using FVM and IFVM

For one dimensional shallow water equations
 we have shown the numerical 
results 
for
1 1
different schemes of FVM. Grid size ∆x = m and time step size ∆t = s
8 3000
have been considered.
The boundary conditions [4] have been taken similar to FDM as

at x = 0 and x = L M = 0.

and η(0, j) = η(0, j) ; η(L − 1, j) = η(L, j)

The initial condition has been taken as

at t = 0, assuming the initial velocity of water is zero, i.e water is at rest position
at t = 0 and therefore M = 0
Similar to the case of FDM we first take the basin depth of water to be zero i.e
the value of h = 0. Another assumption is the maximum and minimum values of η
are 18 m and 20 m respectively.
Thus in case of FVM when the initial velocity of water is zero i.e at t = 0 the
variation of water surface elevation η with different value of x is shown in Fig. 19
using MATLAB.

20
 Figure 19: At t = 0s, graph of η with distance x

The plots for variation of water surface elevation η with different value of
time t in case of upwind interpolation of FVM has been shown in Fig. 20 and
21.

Figure 20: At time t = 0.3333 s Figure 21: At time t = 0.5 s

Similarly the plots for variation of water surface elevation η with different values of
time t in case of central difference interpolation of FVM has been shown in Fig. 22
to 27.

Figure 22: At time t = 0.1 Figure 23: At time t = 0.16667

21
 Figure 24: At time t = 0.3333 Figure 25: At time t = 1

Figure 26: At time t = 1 Figure 27: At time t = 1.6666

As per Eq. (20) the minimum and maximum value of basin depth i.e h =
5 m and h = 95 m respectively have been considered. The behaviour of η with
distance is cited in Fig. 28 and Fig 29.

Figure 28: At t = 0.1s, graph of η withFigure 29: At t = 0.1s, graph of η with

distance x when h = 5m distance x when h = 95m

22
 Figure 30: graph of η with distance x

Numerical results using IFVM

In this case the values of η and M , are in an interval as such the maximum and
minimum values of η are taken in the interval as [17.5 18.5] and [19.5 20.5] respec-
tively. Also, the value of M in the interval has been as [0, 0.5].
After taking the above values in the interval we have investigated the variation
of water surface elevation η with distance x with upper interval, lower interval and
crisp value which is shown in Fig. 30.

6 Comparison between FDM and FVM, conclu-

sion and future directions
A comparative study of numerical results of FDM and FVM for one dimensional
shallow water equations have been presented hare.
After investigating the numerical methods FDM and FVM for solution of one
dimensional shallow water equation we have compared the values for water surface
elevation η with distance x and is shown in the Fig. 31 and 32. Corresponding
comparison tables for both the methods with different time are given in Tables 1
and 2.

23
Table 1: At time t = 0.1 comparision of FDM with FVM for the value of water
surface elevation
S.no FDM FVM
1 20.0000 20.0000
2 19.1820 19.6872
3 20.1911 19.5149
4 19.1784 19.3509
5 20.4563 19.1330
6 19.5764 18.8960
7 20.5701 18.6437
8 20.4156 18.4385
9 19.9754 18.2493
10 18.0000 18.0000

Table 2: At time t = 0.1 for FVM and at time t = .1667 for FDM the value of water
surface elevation
S.no FDM FVM
1 20.0000 20.0000
2 19.5752 19.6872
3 19.7183 19.5149
4 19.4802 19.3509
5 19.3298 19.1330
6 19.2880 18.8960
7 18.7672 18.6437
8 18.8293 18.4385
9 18.1797 18.2493
10 18.0000 18.0000

24
 Figure 32: At time t = 0.1 s for FVM and
Figure 31: At time t = 0.1 s
t = 0.1667s for FDM

Conclusion

Although the results for one dimensional SWEs using finite volume schemes
coincide with the simple schemes of finite difference but advantage of FVM is that
of the meshing scheme. The unique character of finite volume schemes usually
appear when multidimensional problems are solved using unstructured grids.
The most important feature of FVM is that the method is conservative because
the flux entering a given volume is identical to that leaving the adjacent volume and
it is used only when the equations are based on conservation of physical laws. An-
other advantage of the FVM is that it is easily formulated to allow for unstructured
meshes. The method is widely used in computational fluid dynamics (CFD).
Also, in FDM the values of the dependent variables are stored at the nodes only.
In FEM these values are stored at the element nodes. But in FVM, the values of
the dependent variables are stored at the centre of the control volume. In case of
FDM and FEM, conservation of mass, momentum, energy are not ensured at each
cell/control volume. But this is true for FVM. It may be worth mentioning that
FVM takesless time in computation because it converges with less number of con-
trol volumes.

Future Direction

This investigation gives a new idea of the Interval FVM through SWEs and this
can very well be used in future research for better results for other equations obtained
from different applications. The idea may easily be extended to other structured
as well as unstructured problems with various complicating effects. Although this
require more complex forms of interval computation to handle the corresponding
problem.

25
References
[1] Imamura. F, Yalcine. A.C, Tsunami Modeling Manual (Draft) , pp- (1-58),
2006.

[2] George. F, Carrier. T, Harry. Y, Tsunami run-up and drawn-down on a plane

beach. J. Fluid Mech. pp-(79-99), 2005.

[3] Goto. C, Ogawa. Y, Numerical Method of Tsunami Simulation with the Leap-
Frog Scheme.IUGG/IOC Time Project, Manuals and Guides. No.35, Paris, 4
parts, 1997.

[4] Junbo. P, Harvey Mudd College, Numerical simulation Of Wave Propagation

Using the Shallow Water Equations, Web: http://www.math.hmc.edu/ dy-
ong/math164/2007/park/finalreport.pdf [30th April 2014 at 9:31pm], 2007.

[5] Shukla. A, Singh. A. K, Singh. P, A Comparative Study of Finite Volume

Method and Finite Difference Method for Convection-Diffusion Problem. Amer-
ican Journal of Computational and Applied Mathematics ; 1(2): pp- (67-73)
2011.

[6] Boussinesq. F, travaux de M de Saint-Venant,Annales des ponts et chausses

12, pp- (557-595), 1886 .

[7] O. Zikanov. Essential Computational Fluid Dynamics. Wiley India Pvt.Ltd pp-
(86-102), 2012.

[8] H. K. Versteeg, W. Malalasekera. An introduction to computational fluid dy-

namics The finite volume method. Longman Scientific and Technical. pp- (85-
133), 1995.

[9] T. Cebeci, J. R. Shao, F. Kafyeke, E. Laurendeau. Computational Fluid Dy-

namics for Engineers Horizon Publishing. , pp- (141-173), 2005.

[10] C. Dawson, M. Mirabito. Institute for Computational Engineering and Sci-

ences,University of Texas at Austin. 2003.

[11] R. B. Bhat, S. Chakraverty. Numerical Analysis in Engineering, Narosa Pub-

lishing House Pvt. Ltd. 2004.

[12] M. K. Jain, S.R.K. Iyengar, R.K. Jain. Numerical Methods For Scientific and
Engineering Computation, New Age International Limited, Publisher, 2014.

[13] A. batra. Finite Volume Method:basic principles and examples,Department

of Mathematics and Computing, IIT Guwahati, http://www.leb.eei.uni-
erlangen.de/winterakademie/2008/report/content/course02/pdf/0204.pdf,
2008.

26