A nonlinear polar coordinate shallow water model for tsunami computation along North Sumatra and Penang Island

ARTICLE IN PRESS

Continental Shelf Research 27 (2007) 245–257 www.elsevier.com/locate/csr

A nonlinear polar coordinate shallow water model for tsunami computation along North Sumatra and Penang Island Gauranga Deb Roy, Md. Fazlul Karim, Ahmad Izani M. Ismail School of Mathematical Sciences, Universiti Sains Malaysia, 11800 Penang, Malaysia Received 18 February 2006; received in revised form 18 September 2006; accepted 10 October 2006 Available online 8 December 2006

Abstract A nonlinear shallow water model in cylindrical polar coordinate system is developed, using an explicit finite difference scheme with a very fine resolution, to compute different aspects of tsunami at North Sumatra and the adjacent island Simeulue in Indonesia, and the Penang Island in Peninsular Malaysia. The pole of the frame is placed on the mainland of Penang (100.51E) and the model area extends up to the west of Sumatra (87.51E). The model is applied to simulate the propagation of tsunami wave towards North Sumatra, Simeulue and Penang Islands associated with Indonesian tsunami of 26 December 2004. The model is also applied to compute water levels along the coastal belts of those islands. Computed and observed water level data are found to be in good agreement and North Sumatra is found to be vulnerable for very high surges. The computed and observed arrival times of high surges are also in reasonable agreement everywhere. Further studies are carried out to investigate the effect of convective terms and it is found that their effects are insignificant in tsunami propagation and weakly significant for wave amplitude very near to the coast. r 2006 Elsevier Ltd. All rights reserved. Keywords: Tsunami source; Tsunami surge; Tsunami propagation; Sumatra and Penang Islands; Shallow water model; Cylindrical polar coordinates

1. Introduction The Indonesian tsunami of 26 December 2004, having its source near Sumatra, was very severe and caused tremendous loss of life and property along the coastal regions surrounding the Indian Ocean including North Sumatra in Indonesia and Penang Island in Peninsular Malaysia. From this event it is evident that North Sumatra and the Penang Island Corresponding author. Tel.: +60 4 6533964; fax: +60 4 6570910. E-mail addresses: [email protected] (G. Deb Roy), [email protected] (M. Fazlul Karim), [email protected] (A.I.M. Ismail).

0278-4343/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.csr.2006.10.004

are in vulnerable positions for tsunami surges due to a seismic source that may induce static displacement of the sea bottom. It is essential that tsunamis are studied in detail and models be developed on a regional basis for accurate prediction of different aspects of tsunami along the coastal belts. A tsunami is a series of very fast moving long waves generated in a sea by the displacement of water due to undersea landslides, volcanic eruption (Imamura and Imteaz, 1995) and falling of objects from outer space; but mostly, due to seismic fault motions during submarine earthquakes. Tsunamis are characterized as shallow-water waves of extremely long wavelengths; sometimes in excess of 500 km (Yalciner et al., 2005). Although the wave

ARTICLE IN PRESS 246

G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

amplitude is moderate in deep waters (e.g. 0.5–1 m), the tsunami wave slows down and the wave height increases rapidly near the shoreline until it breaks. The wave run-up height along a coastal belt might reach several meters above the mean sea level and cause significant damages. Many works on generation, propagation and runup of tsunami have been done so far. Some of them are Imamura and Gica (1996), Titov and Gonzalez (1997), Imteaz and Imamura (2001), Kowalik et al. (2005), Zahibo et al. (2003), Roy and Ismail (2005), Roy and Ismail (2006). Kowalik et al. (2005) has developed a spherical polar coordinate model for global tsunami computation. This model was applied to simulate the Indonesian tsunami of 26 December 2004 in the World Ocean between 801S and 691N. They used the spherical polar shallow water model incorporating a very fine mesh resolution of 1 min with about 200 million grid points. After the generation of the initial tsunami wave, the disturbance propagates as long waves with a very high speed. So, immediately after the genera-

tion of the tsunami, the people of the surrounding coastal belts must be warned about the arrival time and intensity of the wave amplitude. Estimation of these can best be carried out through a numerical model. Since the time lag between the generation and inundation along a coastal belt is small, any global model, which takes long time for simulation, is not effective for real time simulation. For this purpose a regional model may be more effective and efficient, as it takes less time to simulate and it is possible to incorporate very fine resolution in the numerical scheme. Thus we may have more detailed information about intensity and other aspects of the tsunami around the region of interest. Roy et al. (1999) developed a linear polar coordinate shallow water model to compute tide and surge due to tropical storms along the coast of Bangladesh. Haque et al. (2005) modified the model of Roy et al. (1999), to ensure a resolution that is fine to coarse in the positive radial direction. A logarithmic transformation is used so that in the transformed domain the resolution becomes

Fig. 1. (a) Model domain with the pole of the coordinate system and boundaries. (b) Map of Penang State in Malaysia including Penang mainland and Penang Island.

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

uniform. In the present study, a completely nonlinear shallow water model is developed in cylindrical polar coordinates by setting the pole at the main land of Penang (100.51E) and extending the model region up to west of Sumatra Island (87.51E) (Fig. 1a). The Penang State in Peninsular Malaysia, consisting of the Penang Mainland and the Penang Island, is shown in Fig. 1b. An uneven resolution, fine to coarse in the positive radial direction, is considered and following Haque et al. (2005), a transformation of radial coordinate is applied so that the grid system in the transformed domain becomes uniform. The shallow-water equations and the boundary conditions are transformed into the new domain where an explicit finite difference scheme is used to solve these equations and boundary conditions. This model is used to simulate the propagation of Indonesian tsunami wave towards Sumatra, Simeulue and Penang islands and to estimate the water levels along their coastal belts.

The parameterizations of Fr and Fy are done by conventional quadratic law: qffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffi F r ¼ rC f vr v2r þ v2y and F y ¼ rC f vy v2r þ v2y , (4) where r is the sea water density and Cf is the coefficient of friction. 2.2. Boundary conditions For a closed boundary the normal component of velocity is considered as zero. The radiation type of boundary conditions are used for open boundaries which allow the disturbance created within the analysis area to go out of the area. The analysis area is bounded by the radial lines y ¼ 01, y ¼ Y through O and the circular arc r ¼ R (Fig. 1). Following Roy et al. (1999) the northern, southern, and western open sea boundary conditions are, respectively, given by pffiffiffiffiffiffiffiffi vy þ g=hz ¼ 0 along y ¼ 0, (5)

2. Mathematical formulation 2.1. Vertically integrated shallow water equations To derive the vertically integrated shallow water equations, a system of cylindrical polar coordinates is used in which the pole, O, is in the undisturbed level of the sea surface (the ry-plane) and Oz is directed vertically upwards. The displaced position of the free surface is z ¼ z(r, y, t) and position of the sea floor is z ¼ h(r, y ) so that the total depth of the fluid layer is z+h. The vertically integrated nonlinear shallow water equations are qz 1 q 1 q þ ½rðz þ hÞvr  þ ½ðz þ hÞvy  ¼ 0, qt r qr r qy

(1)

qvr qvr vy qvr qz Fr þ vr þ  fvy ¼ g  , qr rðz þ hÞ qt qr r qy

(2)

qvy qvy vy qvy g qz Fy  þ vr þ þ fvr ¼  , r qy rðz þ hÞ qt qr r qy

(3)

where vr is the radial component of velocity of the sea water, vy the tangential component of velocity of the sea water, Fr the radial component of frictional resistance at the sea bed, Fy the tangential component of frictional resistance at the sea bed, f the Coriolis parameter ¼ 2O sin j, O the angular speed of the earth, j the latitude of the location, and g the acceleration due to gravity.

247

vy 

pffiffiffiffiffiffiffiffi g=hz ¼ 0

along y ¼ Y,

(6)

vr 

pffiffiffiffiffiffiffiffi g=hz ¼ 0

along r ¼ R.

(7)

2.3. Transformation for uneven resolution along radial direction The polar coordinate system automatically ensures finer resolution along tangential direction near the region of the pole. By setting the pole suitably at the location where fine resolution is required, a uniform grid of size Dy is generated in the tangential direction by a set of radial lines through the pole. The arc distance between any two successive radial lines decreases towards the pole and increases away from the pole. Thus, in terms of arc distance, an uneven resolution is achieved in the tangential direction although uniform grid size Dy is used. To achieve uneven resolution along radial direction, fine to coarse in the positive radial direction, according to Haque et al. (2005), the following transformation is used:   r Z ¼ c ln 1 þ , (8) r0 where r0 is a constant of the order of total radial distance and c is a scale factor. From this

ARTICLE IN PRESS 248

G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

transformation we obtain a relationship between Dr and DZ which is as follows: r þ r0 DZ. (9) c This relation shows that, for a constant value of DZ, Dr will increase with increase of r. So, we obtain uneven resolution (fine to coarse) in the radial direction in the physical domain while in computational domain the resolution remains uniform. The Jacobian of the given transformation is as follows:   qr qy    qZ qZ   r0 eZ=c 0  r   c  0 J ¼  qr qy  ¼  (10)  ¼ eZ=c a0  qy qy   0 c 1

Dr ¼

and the operator for the derivative is given by q ceZ=c q  . qr r0 qZ

(11)

Using transformation (8) in Eqs. (1)–(3) we obtain the following transformed equations: o qz ceZ=c q n þ ðz þ hÞðeZ=c  1Þvr Z=c qt r0 ðe  1Þ qZ  1 q  ðz þ hÞvy ¼ 0, þ Z=c r0 ðe  1Þ qy qvr ceZ=c qvr vy qvr þ vr þ  fvy qt r0 qZ r0 ðeZ=c  1Þ qy qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 gceZ=c qz C f vr vr þ vy  , ¼ ðz þ hÞ r0 qZ qvy ceZ=c qvy vy qvy þ vr þ þ fvr Z=c qt r0 qZ r0 ðe  1Þ qy qffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 g qz C f vy vr þ vy  . ¼ ðz þ hÞ r0 ðeZ=c  1Þ qy

ð12Þ

Dy, between any two consecutive straight grid lines through O is constant; whereas the distance between any two consecutive circular grid lines, Dr, increases in the positive radial direction. After transformation (8), both Dy and DZ become uniform. In the transformed domain, the discrete grid points (Zi, yj) are defined by Zi ¼ ði  1ÞDZ;

i ¼ 1; 2; 3; . . . ; M,

(15a)

yj ¼ ðj  1ÞDy;

j ¼ 1; 2; 3; . . . ; N,

(15b)

where M ¼ 878 and N ¼ 307; also Dy ¼ (110/306)1 and DZ ¼ 1/1077 so that in the computational domain y ranges from 01 to Y ¼ 1101 and Z ranges from 0 to 877/1077. Although DZ is a constant, Dr increases with the increase of r according to Eq. (9) and varies from 0.8 to 1.9 km and the total radial distance is RE1160 km. The sequence of time is given by tk ¼ kDt;

k ¼ 1; 2; 3; . . . .

(16)

Although in the physical domain N grid lines meet at the Pole, in the computational domain this point is considered as N distinct grid points (note that it is generated automatically). Since the pole is set at the land, where no computation is done, there is no problem of instability during numerical computation. 3.2. Numerical scheme

ð13Þ

ð14Þ

The boundary conditions (5)–(7) remain unchanged under the transformation. 3. Grid generation and numerical scheme 3.1. Grid generation In the physical domain the grid system is generated through the intersection of a set of straight lines given by y ¼ constant through the pole O (5122.50 N, 1001300 E) and concentric circles, with centre at O, given by r ¼ constant. The angle,

The governing Eqs. (12)–(14) and the boundary conditions (5)–(7) are discretized by finite-difference (forward in time and central in space) and are solved by a conditionally stable semi-implicit method using a staggered grid system, similar to the Arakawa C, as described in Roy et al. (1999). For numerical stability, the velocity components in Eqs. (13) and (14) are discretized in a semi-implicit manner. For example, in the last term of Eq. (13) qffiffiffiffiffiffiffiffiffiffiffiffiffiffi the time discretization of vr v2r þ v2y is implemented qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi as vlþ1 ðv2r þ v2y Þl where superscript l denotes values r of the velocity component at present time step and l+1 denotes the same at the advanced time step. Along the closed boundary, the normal component of velocity is considered as zero, and this is easily achieved through appropriate stair step representation. The time step is taken as 5 s that ensures the CFL stability criterion of the numerical scheme. The values of the friction coefficient is taken as uniform (C f ¼ 0:0033) throughout the physical

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

domain. The depth data used in this study are collected from the Admiralty bathymetric charts.

The generation mechanism of the 26 December 2004 Indonesian tsunami is mainly due to the static sea floor uplift caused by abrupt slip at the India/ Burma plate interface. The extent of the rupture is estimated to be 1200–1400 km from the epicenter along the fault line (Ammon et al., 2005; Hirata et al., 2006; Ohta et al., 2006; Tanioka et al., 2006). A detailed description of the estimation, based on Okada (1985), of the extent of earthquake rapture (between 93–971E and 3–111N) along with the maximum uplift (507 cm) and subsidence (474 cm) of the seabed has been reported in Kowalik et al. (2005). According to this estimation, the tsunami source is elongated along the fault zone from southeast to north-west with uplift to subsidence from west to east. According to the supporting online materials of Ammon et al. (2005) Fig. S2, the extent of the rupture zone is between 92–971E and 2–131N with maximum uplift to subsidence of sea bed within 75 m. According this estimation, the uplift values are 1–5 m across the region with dimension of 100 km  500 km from epicenter near 31N to about 81N. The intensity of uplift/subsidence is least in the region 81N to about 131N. Tanioka et al. (2006) computed the rupture zone between 92–971E and 2–13.51N; the estimated subsidence to uplift are within 3 to 8 m. This estimation is based on tsunami waveform inversion. Based on the above information we consider the source, extended along the fault line between 92–971E and 2–13.51N with maximum rise of 5 m and maximum fall of 4.75 m of the sea surface, as the initial condition in our simulation. Relative to Sumatra and Peninsular Malaysia, the source is in the form of sea level trough to crest. Fig. 2 shows the vertical crosssection of the source along 151st radial grid line. Other than source region the initial sea surface elevations are taken as zero everywhere. The initial radial and tangential components of velocity are also taken as zero throughout the domain. 5. Results and discussions The effects of Indonesian tsunami 2004 along the coasts of North Sumatra, Simeulue and Penang Islands are investigated. Numerical simulation of this potential tsunami is performed in the frame-

5 4 3

elevation in m

4. Initial condition (tsunami source generation)

249

2 1 0 -1 -2 -3 -4 -5

100 200 300 400 500 600 700 800 900 1000 1100

distance in km Fig. 2. Vertical cross-section of the source associated with the Indonesian tsunami 2004 along 151st radial gridline; the distance along the gridline is measured from the pole.

work of nonlinear shallow-water equations in cylindrical polar coordinate system. Wave propagation from the source is computed and the water levels along the coastal belts of Sumatra and Penang Islands are estimated. 5.1. Propagation of tsunami towards Penang and Sumatra islands Fig. 3 shows the contour plot of time, in minutes, for attaining +0.1 m sea level rise at each grid point towards Penang Island. Thus, considering the 0.1 m sea level rise as the arrival of tsunami, it is found that after generation of the initial tsunami wave, the disturbance propagates gradually towards the Penang Island and reaches the west coast at 240 min. Then it gradually proceeds towards the east surrounding the island and reaches the east coast between 280 and 300 min. The propagation time of tsunami along north coast is between 240 and 300 min and the same along the south coast is between 240 and 280 min. The arrival time of tsunami at Batu Ferringi (at north coast) is approximately 270 min. Fig. 4 shows the contour plot of propagation of tsunami towards North Sumatra. The disturbance propagates towards the coastal belts of Sumatra and reaches Banda Aceh (see Fig. 1a for the position) in 20 min. The propagation time along the west coast of North Sumatra is between 20 and 70 min and the same along the east coast is between 20 and 250 min (Fig. 4).

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

250

6°N |

6.5°N 100.5°E

5°N | 28 300 0 2 60

N

4°N 100.5°E

300 280

Penang

100°E -

260

- 100°E

240 220

200

99.5°E -

- 99.5°E

180

160

99°E 6.5°N

|

|

6°N

5°N

99°E 4°N

Fig. 3. Tsunami propagation time in minutes towards Penang (sea level rise of +0.1 m is considered as the arrival of tsunami).

6° N |

7° N 99 ° E

4° N |

N 98 ° E -

110

130 170

210

2° N 99 ° E

250

- 98 ° E

90 70

North Sumatra

50

70

96 ° E 30

Banda Aceh

30 10

10

94 ° E 7° N

- 96 ° E

50

|

|

6° N

4°N

94 ° E 2° N

Fig. 4. Tsunami propagation time in minutes towards North Sumatra (sea level rise of +0.1 m is considered as the arrival of tsunami).

5.2. Computed time series of water levels along Penang and North Sumatra The computed time series of water surface fluctuations at four locations of Penang Island are shown in Fig. 5. The maximum amplitude at Batu Ferringi (at north coast; see Fig. 1b for the position) is approximately 3.3 m (Fig. 5a). At 3 h 40 min after the generation of initial tsunami wave at the source, the water level starts decreasing (withdrawal of water) and reaches the level of 1.2 m. Then the water level increases continuously to reach the level of 1.2 m (first crest) at 4 h 10 min before going down again. The water level oscillates and this oscillation continues for several hours. The computed water

level at Tanjung Tokong (at north coast), shows the same pattern as that of Batu Ferringi with the maximum level of approximately 1.2 m (Fig. 5b). At Pantai Aceh and Pasir Panjang (at west coast) the computed maximum elevations are found to be 3.2 and 3.3 m, respectively (Figs. 5c and d). Fig. 6 depicts the computed time series of water surface fluctuations at four locations on the east and west coasts of North Sumatra and at Simeulue. The rupture area is very near to Banda Aceh and the coseismic vertical displacement induced by rupture extends from the sea bed to the coastline of Banda Aceh. So the coastline of Banda Aceh experiences an instantaneous subsidence of 1.8 m and hence an instantaneous fall of water level of 1.8 m (Fig. 6a).

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

b 5

5

4

4

3

3 elevation in m

elevation in m

a

2 1 0

2 1 0

-1

-1

-2

-2

-3

0

1

2

3

4 5 6 time in hr

7

8

9

-3

10

c

0

1

2

3

4 5 6 time in hr

7

8

9

10

1

2

3

4 5 6 time in hr

7

8

9

10

d 5

5

4

4

3

3 elevation in m

elevation in m

251

2 1 0

2 1 0

-1

-1

-2

-2

-3

0

1

2

3

4 5 6 time in hr

7

8

9

10

-3

0

Fig. 5. Time series of computed elevation at coastal locations of Penang Island associated with the Indonesian tsunami 2004: (a) Batu Ferringhi (b) Tanjung Tokong (c) Pantai Aceh (d) Pasir Panjang.

It is seen in Fig. 6a that, the maximum water level near the coastline of Banda Aceh is 12 m (Fig. 6a); but the actual maximum water level along this coastal belt would be (12 m+1.8 m) ¼ 13.8 m, because of the subsidence of the coastline. The time series of water surface fluctuation at Meulaboh shows that the first crest of amplitude 8.8 m is attained at 45 min after the generation the initial tsunami wave (Fig. 6b). The computed amplitude of the first crest at Medan is approximately 3.1 m and its arrival time is 4 h and 20 min (Fig. 6c). At Simeulue Island the maximum elevation is found to be 6.6 m and the wave arrived after 30 min (Fig. 6d).

5.3. Estimation maximum water level Fig. 7 shows the curves of maximum elevations along four coasts of Penang Island. At the west coast the highest elevation is attained at the southern part (Fig. 7a); as much as 4 m surge is attained at the west coast. At each of the north and south coasts the highest elevation is along the west part (Fig. 7c and d); the maximum surges attained are 4.8 and 3.2 m, respectively. On the other hand, the surge intensity is less at the east coast having 1.4 m as the maximum surge along the south part (Fig. 7b).

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

252

b 14 12 10 8 6 4 2 0 -2 -4 -6

10 8 6 elevation in m

elevation in m

a

2 0 -2 -4 -6

0

1

2

3

4 5 6 time in hr

7

8

9

10

c

0

1

2

3

4 5 6 time in hr

7

8

9

10

0

1

2

3

4 5 6 time in hr

7

8

9

10

d 4 3 2

elevation in m

elevation in m

4

1 0 -1 -2 -3

0

1

2

3

4 5 6 time in hr

7

8

9

10

8 7 6 5 4 3 2 1 0 -1 -2 -3 -4

Fig. 6. Time series of computed elevation at coastal locations of North Sumatra and Simeulue Island associated with the Indonesian tsunami 2004: (a) Banda Aceh (b)Meulaboh (c) Medan (d) Simeulue.

Fig. 8 depicts the similar results along west and east coasts of North Sumatra. The surge intensity is found to be highest at the west coast; the maximum elevation in excess of 22 m (Fig. 8a) is attained at the north of Meulaboh. The maximum elevation along the east coast is found to be 16 m and the elevation at Medan is 6 m (Fig. 8a). 5.4. Arrival time of maximum surge The arrival time of the maximum surge at a given location is an important parameter in the tsunami prediction and warning. The computed time of

attaining +0.1 m above mean sea level (MSL) surrounding the Penang Island and North Sumatra are shown in Figs. 3 and 4, respectively. The times of attaining maximum elevations around these islands are also computed and it is found that this time at each location is 10–15 min later than the time of attaining +0.1 m above MSL. Thus, the time of attaining maximum surge along the west coast of Penang is 250 min (4 h 15 min) and that along the east coast is between 290 and 310 min. The time of arriving maximum surge along the north coast is between 250 and 310 min, whereas the same along the south coast is between 250 and 290 min.

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

b 5

5

4

4 elevation in m

elevation in m

a

3 2 1 0

3 2 1

0

20

40 60 80 100 west coast : north-south

0

120

c

0

20

40 60 80 100 east coast : north-south

0

10

20 30 40 south coast : east-west

120

d 5

5

4

4 elevation in m

elevation in m

253

3 2 1 0

3 2 1

0

20

40 60 80 north coast : east-west

100

0

50

60

Fig. 7. Maximum elevation along the four coasts of Penang Island; the numbers along horizontal axis indicate coastal locations: (a) (west coast) 32—Pantai Aceh, 80—Sungai Burung, 120—Pasir Panjang; (b) (east coast) 20—George Town, 112—Batu Maung; (c) (north coast) 16—Tanjung Tokong, 36—Tanjung Bunga, 76—Batu Ferringi; (d) (south coast) 10—Gertak Sanggul, 50—South of Bayan Lepas.

The maximum surge at Batu Ferringi (at north coast) occurred at approximately 280 min (4 h 40 min). The authors of this paper made a survey to estimate the water levels at different locations of Penang Island and time of arrival of high waves along the coastal belts. On the basis of the authors’ survey and the reports in the newspaper, the arrival time of high surge along the west coast was approximately at 1310 and the same at Batu Ferringi was approximately at 1330 (Malaysian time). The rupture time was 0859 (Malaysian time) and so the propagation time at west coast and at Batu Ferringi were 4 h 11 min and 4 h 31 min, respectively. From the USGS website [http://staff.

aist.go.jp/kenji.satake/Sumatra-E.html (Tsunami travel time in hours for the entire Indian Ocean)] it is found that, the arrival time of maximum surge at Penang Island is 4 h after generation of the tsunami wave. Thus the computed and observed arrival time of high tsunami surge at Penang Island is in reasonable agreement with the observations. The time of attaining the maximum surge along the west coast of North Sumatra is between 30 and 80 min and the same along the east coast is between 30 and 250 min. The computed arrival times of maximum surges at Banda Aceh, Meulaboh, Medan and Simeulue are, respectively, 30, 45, 250 and 30 min. According to the survey report the tsunami

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

254

b 24

24

20

20 elevation in m

elevation in m

a

16 12 8 4 0

16 12 8 4

0

50

100 150 200 west coast : north-south

250

0

0

50

100 150 200 250 east coast : north-south

300

Fig. 8. Maximum elevation along the two coasts of North Sumatra; the numbers along horizontal axis indicate coastal locations: (a) (west coast) 15—West Banda Aceh, 50—North Meulaboh, 70—Meulaboh, 90—Meulaboh Aronghan; (b) (east coast) 40—North of Banda Aceh, 80—Singli, 200—Langsa, 270—Medan.

waves have arrived near Banda Aceh in half an hour and have reached at Meulaboh, Medan and Simeulue at 40, 240 and 30 min, respectively (Yalciner et al., 2005). Thus there is a matching between the model results and the data of observations. 5.5. Comparison between model results and observed data The computed water levels along the North Sumatra and Penang Island will now be compared with the observed water level data. A post-tsunami survey has been carried out by the authors of this paper at the coastal regions of Penang Island to estimate the maximum water level above the MSL associated with Indonesian tsunami 2004. A Turkish–Indonesian–USA team (Yalciner et al., 2005) made another survey along the coasts of North Sumatra to estimate the water levels and runup heights for the same event. Observed and computed water levels at different locations of Penang Island (Table 1) are in a reasonable agreement except at Tanjung Tokong and Tanjong Bunga, where the computed water level is less than the observed value. The computed surge levels at the west part of the north coast appears to be higher (Fig. 7c) than the observed values; the reason for which should be investigated. On the other hand, the computed results confirm the

Table 1 Observed and computed water levels above the mean sea level at different locations of Penang Island Location name

Observed water level (m)

Batu Ferringhi (north penang) Taluk Bahang Tanjung Tokong Tanjong Bunga Pasir Panjang Kuala Pulau Betung

3.3 42.0 2.8 3.0 42.0 2.0

Computed water level (m) 3.3 3.76 1.89 1.47 3.3 2.49

Source: The authors conducted this survey.

Table 2a Observed and computed water levels above the mean sea level at different locations of North Sumatra Location name

Observed maximum tsunami amplitude near the shore (m)

Simeulue Island Meulaboh Medan Banda Aceh Meulaboh Aronghan St 176, WPT 17

43

Computed water level (m)

6.6

1.5 42.5 — 15

8.8 3.1 13.8 16.94

415

19.10

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257 Table 2b Observed and computed arrival time of tsunami surge at different locations of North Sumatra and Simeulue Islands Location name

Observed arrival time of tsunami amplitude near the shore (min)

Computed arrival time of maximum surge (min)

Simeulue Island Meulaboh Medan Banda Aceh Meulaboh Aronghan St 176, WPT 17

E30 E40 E240 E40

30 45 250 30 35

E40

33

Studies are carried out to investigate the effect of nonlinear terms of the momentum equations. Simulations are performed with the corresponding

b 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2

l60 nl60

elevation in m

elevation in m

observed directivity of tsunami amplitude at different locations of North Sumatra (Table 2a). The water levels at different locations are found to be in good agreement with the data of observations. The computed and observed times of attaining maximum surges at some coastal locations along the North Sumatra are also shown in Table 2b. The computed and observed times at those locations are almost identical. 6. Investigation on the effect of nonlinear terms

a

200

400 600 800 1000 distance in km from the pole

c

3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2

ll20 nll20

200

400 600 800 1000 distance in km from the pole

d 3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2

l220 nl220

elevation in m

elevation in m

255

200

400 600 800 1000 distance in km from the pole

3.5 3 2.5 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2

l264 nl264

200

400 600 800 1000 distance in km from the pole

Fig. 9. Computed disturbance patterns associated with linear and non-linear models along 151st radial gridline (the distance along the gridline is measured from the pole): (a) at 60 min, (b) at 120 min, (c) at 220 min, and (d) at 264 min (time of hitting the coast).

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

256

linear model by excluding the convective terms and the results are compared with those of the nonlinear model. Fig. 9 shows the propagation of the disturbance along the 151st radial grid line starting from the western open boundary to the west coast of Penang Island (1160 to 60 km from the pole) for both the models. Solid and dotted lines represent the disturbance patterns computed through nonlinear and linear models, respectively. This figure shows the disturbance pattern at 60, 120, 220 and 264 min (time of attaining maximum water level at the coast) along the 151st radial grid line. Considering the time first sea surface rise as the time of arrival of tsunami, it is seen that tsunami surge reaches the locations 335, 210 and 75 km from the pole at 60, 120 and 220 min, respectively for both the cases. Thus, the computed propagation pattern associated with linear and nonlinear models are almost identical. This indicates the insignificance of convective terms and confirms that their effect is negligible in computing tsunami travel time. However, the maximum water levels at 264 min at the coast are 3.2 and 3.35 m associated with the nonlinear and linear models, respectively (Fig. 9d). Thus, nonlinear terms are found to be weakly significant for the wave height near the coast. To investigate effect of nonlinear terms in the very shallow region, the portion of the Fig. 9d between 60 and 110 km (that is 50 km from the coast) is zoomed and is shown in Fig. 10a. In this region the depth of the sea varies from 5 to 50 m. This figure

a

7. Conclusions In this study we have developed a nonlinear cylindrical polar coordinates model for the region covering Penang and North Sumatra and is used to compute different aspects of Indonesian tsunami 2004. Based on the discussions above, the following conclusions can be made: The model is simulating the propagation of tsunami wave from the source with reasonable accuracy. The computed water levels along the coastal belts of the Islands are in reasonable agreement with data of observations, except at some locations of Penang Island. Factors such as varying bathymetry and coastline geometry can affect the amplitude of the waves along the coastal belts. The computed arrival times of tsunami at

b

3.5 3

shows that, the surge amplitude increases very fast near the coast having the depth between 5 and 20 m. This figure also demonstrates that, the water level at the coast is 3.2 m (solid line) for the nonlinear model, whereas the same is 3.35 m (dotted line) for the linear model. The wave shapes are almost identical for both the models within 20 to 50 m depth. Fig. 10b shows the peak elevations along the 151st radial grid line between 60 and 110 km (50 km from the coast) associated with both the models. The variation in the elevations between the results of two models is observed up to 30 km from the coast and at the coast the elevations associated with nonlinear and linear models are 3.2 and 3.5 m.

l264 nl264

3 elevation in m

elevation in m

2.5 2 1.5 1 0.5 0

l151 nl151

2.5 2 1.5 1 0.5

-0.5 -1 60

3.5

70 80 90 100 distance in km from the pole

110

0 60

70 80 90 100 distance in km from the pole

110

Fig. 10. Comparison of two phenomena associated with non-linear (solid line) and linear models(dotted line) within 50 km from the west coast along 151st radial gridline: (a) disturbance patterns at 264 min, and (b) maximum water levels.

ARTICLE IN PRESS G. Deb Roy et al. / Continental Shelf Research 27 (2007) 245–257

different locations of those Islands are in reasonable agreement with observed arrival times. Nonlinear terms are found to be insignificant for wave propagation (travel time); however they are weakly significant for the wave amplitude near the coast. So the nonlinear model should be applied to compute wave amplitude near the coast. This model may be applied to compute the tsunami phenomena to any localized coastal region where a tsunami hazards exist. Acknowledgements This research is supported by a short-term grant of Universiti Sains Malaysia (USM) and the authors acknowledge the support of the USM short-term grant. References Ammon, J.C., Ji, C., Thio, H., Robinson, D., Ni, S., Hjorleifsdottir, V., Kanamori, H., Lay, T., Das, S., Helmberger, D., Ichinose, G., Polet, J., Wald, D., 2005. Rupture process of the 2004 Sumatra–Andaman Earthquake. Science 308, 1133–1139. Haque, M.N., Miah, S., Roy, G.D., 2005. Polar coordinate shallow water model with radial stretching for the coast of Bangladesh. SUST Studies (J. Shah Jalal Univ. Sc. & Tech., Sylhet, Bangladesh) 6, 1–13. Hirata, K., Satake, K., Tanioka, Y., Kuragano, T., Hasegawa, Y., Hayashi, Y., Hamada, N., 2006. The 2004 Indian Ocean tsunami: tsunami source model from satellite altimetry. Earth Planets Space 58, 195–201. Imamura, F., Gica, E.C., 1996. Numerical model for wave generation due to subaqueous landslide along a coast—A case of the 1992 Flores tsunami, Indonesia. Science of Tsunami Hazards 14 (1), 13–28. Imamura, F., Imteaz, M.A., 1995. Long waves in two layers: governing equations and numerical model. Science of Tsunami Hazards 13 (1), 3–24.

257

Imteaz, M.A., Imamura, F., 2001. A non-linear numerical Model for stratified tsunami waves and its application. Science of Tsunami Hazards 19 (3), 150–159. Kowalik, Z., Knight, W., Logan, T., Whitmore, P., 2005. Numerical modeling of the global tsunami: Indonesian tsunami of 26 December 2004. Science of Tsunami Hazards 23 (1), 40–56. Ohta, Y., Meilano, I., Sagiya, T., Kimata, F., Hirahara, K., 2006. Large surface wave of the 2004 Sumatra–Andaman earthquake captured by the very long baseline kinematic analysis of 1-Hz GPS data. Earth Planets Space 58, 153–157. Okada, Y., 1985. Surface deformation due to hear and tensile faults in a half-space. Bulletin of the Seismological Society of America 75 (4), 1135–1154. Roy, G.D., Ismail, A.I.M., 2005. An investigation on the propagation of 26 December 2004 tsunami waves towards the west coast of Malaysia and Thailand using a Cartesian coordinates shallow water model. Proceedings of the International Conference in Mathematics and Applications. Mahidol University, Thailand, pp. 389–410. Roy, G.D., Ismail, A.I.M., 2006. Numerical modeling of tsunami along the coastal belt of Penang using a polar coordinate shallow water model. Far East Journal Applied Mathmatics 23 (3), 241–264. Roy, G.D., Humayun, K.A.B., Mandal, M.M., Haque, M.Z., 1999. Polar coordinate shallow water storm surge model for the coast of Bangladesh. Dynamics of Atmospheres and Oceans 29, 397–413. Tanioka, Y., Yudhicara, Kususose, T., Kathiroli, S., Nishimura, Y., Iwasaki, S.I., Satake, K., 2006. Rupture process of the 2004 great Sumatra–Andaman earthquake estimated from tsunami waveforms. Earth Planets Space 58, 203–209. Titov, V.V., Gonzalez, F.I., 1997. Implementation and Testing of the Method of Spliting Tsunami (MOST) Model. NOAA Technical Memorandum ERL, PMEL -112, Contribution No. 1927 from NOAA/Pacific Marine Environment Laboratory, p11. Yalciner, A.C., Karakus, H., Ozer, C., Ozyurt, G., 2005. Understanding the generation, propagation, near and far field impacts of tsunamis. Lectures Notes, Middle East Technical University, Ankara, Turkey. p135. Zahibo, N., Pelinovsky, E., Kurkin, A., Kozelkov, A., 2003. Estimation of far-field tsunami potential for the Caribbean coast based on numerical simulation. Science of Tsunami Hazards 21 (4), 202–222.