STUDY ON APPROPRIATE MODELING OF TSUNAMISIN MALAYSIA FOR RISK EVALUATION
Zaty Aktar binti Mokhtar*

Supervisor:

MEE06025

Fumihiko Imamura**
Shunichi Koshimura**

ABSTRACT
In order to design a tsunami warning system in a region, there are three major issues in developing thedata base or estimating tsunami risk; the first is the assumption of tsunami sources by selecting thepotential and past earthquakes and others, the second is accuracy and stability on the numerical modelwith the computational regions for tsunami propagations, the third is output such as water level andvelocities from the model to provide the warning and to estimate the risk. In this study, a numericalsimulation of tsunamis in the Sabah coast is conducted. Two types of models are developed in this study:(1) a numerical model of the generation and trans-oceanic propagation of tsunamis using linear theory inspherical coordinate system and (2) a numerical model based on the nested grid system in Cartesian
coordinate system using linear and non-linear theory with a two different spatial grid sizes. Six cases offault ruptures are considered in the Manila trench for the earthquakes with magnitudes of 9.0, 8.5 and 8.0
which are corresponding to the earthquakes with the return periods of 667, 205 and 63 years, respectively.
In the nested grid system model, the linear shallow water wave theory in spherical coordinate system isused for tsunami simulation in the large area covering Southeast Asia while the non-linear shallow water
wave theory in Cartesian coordinate system is used for tsunami simulation in the Sabah coast region. It isfound that the tsunamis arrives the Northern part of Sabah in approximately 2 hours after an earthquake
occurred and the maximum tsunami height calculated is more than 1 m for the Mw 9.0 earthquake. Finally,a few computational instability problems are introduced and discussed in this paper.
Keywords: Numerical Model, Tsunami Simulation, Nested Grid SystemINTRODUCTION
The Sumatran mega-thrust earthquake on 26 December 2004 has awakened specially the affected
countries for the need of an early warning system and re-evaluation of the tsunami hazard map in theregion. Malaysia is located in the area which can be affected by the tsunamis generated by the earthquakes
in the seismically-active western part of the Philippines. Therefore, further research should be emphasized
for this possibility of the future tsunami events in the region. In this study, the simulation of tsunamis in
the Sabah coast area using nested grid system method is conducted to investigate the tsunami arrival times,
and wave heights at several selected points. The numerical simulation is conducted using TUNAMI codes(Imamura et al., 2006) with modifications to simultaneously solve tsunami propagation in the two regionswith using the Leap-Frog Scheme (Goto et al., 1997).

________________________________* Malaysian Meteorological Department, ** Disaster Control Research Centre, Tohoku University Japan1

EARTHQUAKES IN THE PHILIPPINES AND FAULT PARAMETERS
The geometry of subduction slabs in northern Luzon, Philippines has been studied by Bautista et al.,(2001), where a new model of the subducting slabs of the Eurasian plate beneath the Manila Trench wasproposed. They used the hypocentral and focal mechanism data in order to characterize the tectonic
configuration in the Northern Luzon region. Papazachos et al., (2004) have proposed the empiricalformulas for estimating the subsurface fault length, L, fault width, W and displacement, u from themoment magnitude, Mw. The equations were developed from global earthquakes classified into strike-slip
faults, dip-slip continental faults, and dip-slip faults in subduction regions. Therefore, based on theempirical equations, the static fault parameters for each magnitude which are the length, width, anddisplacement are predicted and summarized in Table 1 below.
Table 1 Predicted fault parameters estimated using empirical eqs. (Papazachos et al., 2004)
Magnitude
(Mw)
8.0
8.5
9.0

Subsurface
rupture length, (L)(km)

Subsurface rupture
width (W)(km)

Displacement,(u)
(m)

162.2

70.8

2.19

305.5

101.2

4.57

575.4

144.5

9.55

In this model, for the dip angle of a fault plane is assumed equal to the one proposed by Bautista et al.
(2001); Dip = 30o and (Strike, Slip) = (0o, 90o).
As the initial condition in this model, I assumed that the vertical displacement of water surface isinstantaneous as the same as vertical displacement of the sea bottom which is calculated by using
Mansinha and Smylie (1971) theory. Six cases of fault ruptures are considered for the earthquakes withmagnitudes of 8.0, 8.5 and 9.0. The return periods of earthquakes with the magnitude Mw 9.0, 8.5 and 8.0
are 667, 205, 63 years, respectively (Ruangrassamee, 2007). Top depth of the fault is assumed 25 km. The
fault planes for all six cases are located on the subducting slab in the Manila Trench as shown in Figure 1.

Mw = 8.0
Mw = 8.5Mw = 9.0
Figure 1 Location of fault plane models for six cases in the Manila Trench area (Sources:Philippine Institute of Volcanology and Seismology, http://www.phivolcs.dost.gov.ph/).
)2

ANALYTICAL SCHEME AND CONDITIONS
Figure 2 depicts the area of computation considered in the analysis. In this paper, the nesting grid system
method is introduced where there are two different grid size regions which are considered in this model. In
order to calculate the trans-oceanic tsunami propagation in the long distance (region 1) , the linear longwave equations with Coriolis force in sphericalcoordinate system are used in this model which isnumerically solved by using finite difference
method. Meanwhile, the domain of computation
for region 2 covering the Sabah coast area is in
Cartesian coordinate system which used the nonlinear shallow water wave equations in numericalcomputation. The sea depth off the shore of thePhilippines is about 4 km to 5 km, while the depth
around Sabah coast area is less than 2 km and quite
shallow. The digital global bathymetry grid data
are extracted from (British Oceanographic DataCentre, 1997) for both regions with spatial grid size
1-arc-minute (about 1850 m) for trans-oceanic
tsunamis calculation (region 1). As for region 2, I
used modified GEBCO data with a resolution of 20
sec (about 616.67 m).
Figure 2 Computational areas in numerical model. In order to stabilize the numerical computation, thetemporal grid sizes (dt) were set to 3.0 sec forregion 1 and 1.0 sec for region 2, which satisfied the stability condition in the analysis. The number ofgrid points are 1381× 1201and 772 × 631 for region 1 and 2, respectively. The computation time is about12 hours.

ANALYTICAL RESULTS AND DISCUSSION
The arrival times and tsunami
heights at several selected points
in the case study area areinvestigated for the tsunami
warning and risk evaluation.
Results of both methods arecompared and it is found that thetsunami amplitudes obtained forthe linear theory model werelarger than the non-linear theorymodel at the output points. It
means that the predictions by thelinear model are overestimated but
it is safer evaluation. This
happened by neglecting bottom
friction and convection terms in
the calculation of tsunami

Sabah

Figure 3 Snapshot of tsunami propagation after 3 hours
using the nested grid system method.

3propagation. Figure 3 shows the snapshot of tsunami wave 3 hours after the earthquake occurred with thetime histories of the water level computed for each output point in the region using nested grid system
method. It is noticed that, the tsunamis arrived at thenorthern part of Sabah approximately in 2 hours after
Vietnaman earthquake Mw 9.0 occurred with the maximum
tsunami height of more than 1 m. In this model, it isobvious that the fault rupture occurred off the shore of
the Philippines for the Mw 9.0 has generated a strong
directivity of tsunami energy to propagate toward
Vietnam coast area and a very small affect to thecoastline of Sabah (see Figure 4). At the western coastof Sabah, the calculated maximum tsunami height isMalaysia
about 1.25 m at a sea depth of 5 m in Labuan output
point. However, it is also noticed that the tsunami
height at some other places are likely more than 2 m.
Therefore, the characteristic of sea bottom topography
is really important and needs to be emphasized, in fact
it is one of the major factors determining the influenceand severity of the tsunami. The effect of earthquake
magnitudes on the tsunami height and arrival time at
Figure 4 Directivity of the tsunami front
several output points are presented in Table 2.
towards Vietnam (linear model).
Table 2 Maximum tsunami heights and tsunami arrival times at selected output points
Mw
9.0
8.5

8.0

Case

Max. Tsunami Height (m)Tsunami Arrival Time (min)Kdt

KKLbnSdkn

Kdt

KKLbn

Sdkn

1

1.15

0.79

1.25

0.16

120

156

170

165

2

0.34

0.19

0.39

0.05

141

172

184

336

30.51

0.19

0.51

0.08

126

158

170

301

40.16

0.05

0.10-157

187

201-

5

0.11

0.06

0.13

0.03

143

174

187

511

6

0.18

0.07

0.12

0.03

130

162

175

458

Kdt: Kudat, KK: Kota Kinabalu, Lbn: Labuan, Sdkn: Sandakan
COMPUTATIONAL PROBLEMS
There are some computational problems occurred during this study effort and some existing
approaches can be used to solve the problems. One example of major instability occurred in thecomputation as shown in Figure 5. We can notice that three sources of instability at the bottom rightcorner of the boundary are indicated with the red circles (see Figure 5(a)). Since the depth at these threesources is estimated around 1 m ~ 2 m and the equation used to calculate the water level at this boundary4

is a linear equation theory, thus, it was impossible to calculate the water level at the boundary in this case.
A large dislocation between depth in region 1 (linear theory –the minimum depth is set up to 10 m only)
and region 2 (non-linear theory – 1 m ~ 2 m minimum depth) yields unrealistic results of water level at theboundary region. This is the main cause of the instabilities and demolishes the whole computation asshown in Figure 5(b).
There are many ways to eliminate this kind of instability, for example I have created a bit wider
new bathymetry data. However, I have to make sure that the depth around the boundary of region 2 ismore than 10 m in order to make appropriate the linear equation used in the boundary of region 1 and toavoid the same problem of instability.

(b)

(a)

Sabah

Figure 5 Initiation of an instability at a bottom right corner of the boundary (leftpanel) which demolished the whole computation (right panel).
CONCLUSIONSIn this study, both linear theory and non-linear theory model based on nested grid system computations oftsunami propagation are conducted in the case study area. The simulations of tsunamis are conducted
appropriate with some conditions, where the results of both methods are compared. It is noticed that thefault ruptures occurred off the shore of the Philippines will generate a strong directivity of tsunami energy
to propagate toward Vietnam coast area and a very small affect to the coastline of Sabah.
The tsunami arrival time and maximum wave amplitude for case 1 to 6 are analyzed at eachselected output points. It is noticed that, the tsunamis arrived earlier at each output point in case 1
compared with another cases 2 to 6. The calculated maximum tsunami heights is 1.25 m in Labuan forcase 1, 0.39 m in Labuan for case 2, 0.51 m in Kudat for case 3, 0.16 m in Kudat for case 4, 0.13 m inLabuan for case 5, and 0.18 m in Kudat for case 6.
However, there are some important points that need to be taken into account after conducted thesesimulation. In order to get much more consistent and better results for tsunami computation, the use ofdetailed and accurate bathymetry data is essential with a small grid size and is much more effective in thenumerical computation.
In this paper, I have examined the computational instability problem that occurred during mystudy effort. A better knowledge is required to modify the existing program which can minimize the error
and to avoid instability in the computation.
On the whole, the tsunami heights are larger in the western part of Sabah than the eastern part.
The lifeline facilities and the residents living near the coastal zone should be well prepared for the disaster.

5

Finally, the structural counter measures using natural protection method (soft structure) such as
mangrove plantation or green belt along the coasts is a good protection in order to reduce the tsunami
threat in the case study area.

ACKNOWLEDGEMENTS
I wish to express my special gratitude to all of Graduate School of Engineering, Tohuku University andProfessor Nobuo Shuto (Research Adviser of ARISH) of Nihon University for their helpful discussionsand valuable comments on the manuscript as well as their guidance during this study. I wish to thankmembers of DCRC, particularly Mr. Hideaki Yanagisawa and Mr. Ikuo Abe for their numerous helps andassistances during the study effort. I am grateful to acknowledge Dr. Bunichiro Shibazaki (Tsunami
Course Leader) and Dr. Yushiro Fujii for their guidance, suggestions and encouragements during mystudy in Japan. I thank all the staff and lecturers in IISEE for providing me very valuable knowledge andmaterials that are really useful for my future needs.
REFERENCES
Bautista, B. C., et al., 2001, Tectonophysics, 279-310.
British Oceanographic Data Centre, 1997, “The Centenary Edition of the GEBCO Digital Atlas (CDROM)”

Goto, C., et al., 1997, IOC Manual, UNESCO, No.35.
Imamura, F., et al., 2006, Tsunami modeling Manual (TUNAMI model), 72p.
Mansinha, L. and D. E., Smylie, 1971, Bulletin Seismological Society America, 61, 1433-1440.
Papazachos B. C., et al., 2004, Bulletin of the Geological Society of Greece, 36, 1482-1489.
Ruangrassamee, A. and N., Saelem, 2007, Fourth International Conference on Urban Earthquake
Engineering, 447-452.

6