Understanding the Simple Model of Tsunami Propagation by SiTProS - - PDF document

Understanding the Simple Model of Tsunami Propagation by SiTProS Model

Wattana Kanbua1*, Somporn Chuai-Aree2

1 Marine Meteorological Center, Thai Meteorological Department, Bagkok 10260, Thailand 2 Faculty of Science and Technology, Prince of Songkla University, Pattani 94000, Thailand

E-mail: watt_kan@hotmail.com* ABSTRACT

In this paper, we use the approximation of shallow water waves to understand the behaviour of a tsunami in a variable depth. We deduce the shallow water wave equation and the continuity equation that must be satisfied when a wave encounters a discontinuity in the sea

• depth. Our procedure also includes a new mathematical model for tsunami generation, propagation, real-time simulation and visualization.

The model is so called SiTProS (Siam Tsunami Propagation Simulator). The model can run for any given regional or global grid with a prescribed topographic dataset as ETOPO2. The grid resolution can be arbitrary in space and time. The propagation can be done on a latitudinal – longitudinal or on a cartesian grid. The finite differences method is used to solve the equation. A tsunami struck the shores of southern Thailand and along the Andaman coast on the December 26, 2004. The hardest-hit area of affected provinces is given based on the shoaling, refraction, diffraction and reflection phenomenon. The SiTProS can run for any given regional or global grid with a prescribed topographic dataset as Etopo2. The finite difference method is used to solve the equation. The SiTProS provides five different computing regions, which are Asia, Europe, Arab (Iran-Kenya), Africa, and Andaman. The SiTProS is designed for interactive simulation and user friendly. This paper has introduced the tsunami propagation simulator model which calculated trends to check the model could predict the tsunami arrival times. We have already compared the result of model with observation time in the case of the December 26, 2004. In general, it was found that the model can predict quite well the timing. The SiTProS as tsunami propagation simulator model has therefore been developed by ourselves, using a fast algorithm for a quick estimate of the tsunami front propagation. We have applied this model to the 26th December event. It shows a very satisfied prediction of the time of the main events. The software is available at http://www.schuai.net/SiTProS.

• 1. INTRODUCTION

On December 26, 2004 at 07:59 am (UTC 00:59 am, JST 09:59 am), a giant earthquake occurred off the west coast of northern Sumatra,

• Indonesia. Its epicenter is shown in Figure 1. Figure 2 shows that the seismic activity in this region is very high as the Pacific Rim.

Its magnitude was reported by some institutes as shown in Table 1. The West Coast Alaska Tsunami Warning Center and the Pacific Tsunami Warning Center issued magnitude 8.0 warnings within fifteen minutes of the earthquake. The magnitude 9.0 which is well- known nowadays was reported nineteen hours later. These revised magnitudes, however, do not mean the centers made mistakes, but indicate the difficulty in analyzing such a giant earthquake in such a short time, even for leading seismologists. This was the fourth largest earthquake in the world since 1900, see Figure 3 and Table 2. Another earthquake which occurred in this region three months later was the seventh largest event. The giant earthquake generated a huge tsunami which was the third largest since 1900, as shown in Table 3. This tsunami hit many countries in the Indian Ocean. With the exception of Indonesia, the Andaman Islands and Nicobar Islands, the tsunami, not the earthquake, caused all of the extensive damage (see Figure 4.). This was the greatest tsunami disaster in history. The tsunami hit the southeast coast of Thailand, which was about 500 km from the epicenter. Because the area has world famous resorts like Phuket Island and the tsunami hit the coast at around high tide (Figure 5), there was a dreadful tragedy. 5,400 people were killed and 3,100 people reported missing due to the tsunami in Thailand. To study this disaster, a field survey was carried out from December 30, 2004 to January 3, 2005 along the southeast coast of Thailand. Further, a numerical simulation was conducted to investigate the source mechanism of the tsunami. In this chapter, those results are reported.

Figure 1. The epicenter of the 2004 Sumatra Earthquake. Figure 2. Seismic activity in the world from 1978 to 2000 [4]. Table 1. A history of reported magnitude of the earthquake [3].

*1 "M" means that the type of magnitude was not shown in the e-mail. *2 Where no issued time was shown in the e-mail, the posted time informed by the institute's SMTP server is used. *3 The origin time of the earthquake is assumed to be 12/26/2004 00:59 UTC by USGS.

Table 2. Largest earthquakes in the world since 1900 [3] Table 3. Largest tsunamis in the world since 1900 [1]

Figure 3. Bathymetry. Figure 4. The tidal change and the tsunami arrival time at Phuket [5].

• 2. FIELD SURVEY AND DATA FOR THE EVENT

The general objective of the field survey on the disaster is to determine the damage caused and to study the factors involved. The magnitude of the damage is a result of a balance between factors on the side of the disaster and factors on the human side. In a tsunami disaster, the former factors imply tsunami height, velocity, hydraulic power, etc. and the latter factors imply preparedness, countermeasures, education, evacuation, etc. Many field surveys are desired to investigate factors on both sides.

Because of the difference in tide levels between when the tsunami arrived and when our measurements were made, the measured tsunami heights were corrected by the method shown in Figure 5. The tsunami arrival times were assumed to be a uniform 10:00 am local time.

Figure 5. The method of tide level correction [5]. In order to evaluate the various calculations an analysis of the arrival times of the waves on the various locations has been performed. Three sources of information have been used such as news reports, radar images from the NASA/NOA satellites and surveys data.

• 3. THE SITPROS MODEL

The SiTProS model stands for “Siam Tsunami Propagation Simulator” model which is Tsunami Propagation Model. The SiTProS has simulated and animated tsunami generation and propagation in a given arbitrary shaped bathymetry. The module for wave equation, describe each variable and parameter. We calculate from tsunami behavior, shallow water equation by defining wave propagation speed.

The model for Tsunami propagation is based on wave equation. The code is programmed by Pascal Programming language. This software is designed for fast computing in Real-Time Simulation and Visualization in 2D domain based on graphical user friendly interface. This model and all simulated results cannot be used for any commercial purposes, but training, warning systems and education are available.

• 4. NUMERICAL SIMULATION

4.1. Model Equation

The numerical code is based on wave equation in Grid computing, we have rectangle grid, 2D grid and coordinate system in Figure 6(a), (b). We calculate from tsunami behavior, shallow water equation by defining wave propagation speed. The wave equation can be written as

2 2 2 2 2 2 2

U U U a t x y ⎛ ⎞ ∂ ∂ ∂ = + ⎜ ⎟ ∂ ∂ ∂ ⎝ ⎠

(1)

where U is wave height, x is spatial grid in x-direction, y is spatial grid in y-direction, a is wave propagation speed, d is water depth.

The shallow water equation can be written as

d g a × =

(2) (a) (b) Figure 6. 2D grid system.

Approximating all second derivatives is using central differences for finite different method (FDM).

2 1 , , 1 , 2 2

) ( 2 t U U U t U

n j i n j i n j i

Δ + − = ∂ ∂

+ − 2 , 1 , , 1 2 2

) ( 2 x U U U x U

n j i n j i n j i

Δ + − = ∂ ∂

+ −

(3)

2 1 , , 1 , 2 2

) ( 2 y U U U y U

n j i n j i n j i

Δ + − = ∂ ∂

+ −

From (3),

1 + n

U

can be calculated as :

, 1 2 , 1, 1, , 1 1 , 1 , , ,

( ) ( 4 ) 2

i j n n n n i j i j i j i j n n n n i j i j i j i j

a t U U U U h U U U U

+ − + − − +

Δ = + + + Δ − + −

(4) where

j i

U , is wave height at point j i,

, i is grid index in x-direction, j is grid index in y-direction, a i,j is wave propagation speed at position i,j, ∆t is time step, and ∆h is grid resolution. From (4), we modified the tsunami propagation speed at the position i,j to their four neighboring positions. It is written as equation (5).

1 2 2 2 , 1, 1, 1, 1, 2 2 2 , 1 , 1 , 1 , 1 , , 1 , ,

( ) ( 4 ) 2

n n n i j i j i j i j i j n n n i j i j i j i j i j i j n n i j i j

t U a U a U h a U a U a U U U

+ − − + + − − + + −

Δ = + + Δ + − + −

(5) The wave propagation term a is calculated from the shallow wave property from (2). Since we use the FDM method to approximate the wave height U, we have modified the wave propagation term in the following form:

, 2 , max i j i j

d a d =

(6) where ai,j is tsunami propagation speed at position i,j, di,j is the water depth at position i,j, and dmax is the maximum water depth in the domain. The time step ∆t is calculated from (7). Its value is depended on the initial tsunami propagation speed V0 (m/s).

3700 ( ) 60 3.6 V t Δ =

(7) The grid resolution ∆h is equal to 3700 meters for Etopo2 and 1850 meters for Etopo1 (∆h=∆x=∆y). The 2 minutes digital bathymetry data published Smith and Sandwell (1997) is used to compute the model. We re-sampled this data and produced divided grid bathymetry data for numerical simulation. The model simulates and visualizes Tsunami

(a) Asia (b) Europe (c) Arab (Iran-Kenya) (d) Andaman Sea (e) Africa

Figure 7. The digital bathymetry for computing grid (a) Asia [50N-20S, 80E-150E], (b) Europe [80N-10N, 25W-45E], (c) Arab (Iran-Kenya) [30N-40S, 30E-100E], (d) Andaman Sea [30N-20S, 60E-110E], and (e) Africa [20N-50S, 40W-30E].

generation and propagation in a given computing grid in five different regions such as Asia, Europe, Arab (Iran-Kenya), Africa, and Andaman on Etopo2.

4.2. Tsunami Source Model

For the present numerical simulations, the initial fault displacement was inferred from one of the two Harvard CMT solutions (Table 4). Two fault models assumed to cover distribution of aftershocks, and these sizes are similar, 500km x

• 200km. The resulting sea floor deformation was computed by the theoretical method [2] (Figure 8). Then crustal rigidity is

11

10 4

−

×

dyne/cm2. In the numerical computation, dynamic fault parameters are considered. That is, sea bottom deformation starts at the point of main shock, and radiates at rates of rupture velocity of 2.5km/s. Rise time at each point is 0.1 of a second. Time for crustal deformation is about 400 seconds, because the length of the whole fault area is larger than 1000km.

Figure 8. Fault parameters, sea bottom deformation, gap among contour lines is 0.5m. [5]

4.3. Algorithm and Fast Computation

We develop and improve the algorithm for fast tsunami propagating calculating based on the moving boundary of computing grid. It means at the beginning of initial fault, the small computing region is computed and the rectangle boundary will be moved and expanded for the next time step. In order to guarantee that the solution is controlled, the progressing wave is bounded by the algorithm. This method can help the calculation to reduce the time consuming. The performance of computing time is acceptable and fast enough for the warning

• system. The following steps are algorithm in this calculation:
• 1. Load computing map by user selection
• 2. Select computing region from a whole map and compute (6), (7)
• 3. Setup the initial tsunami parameter effected by fault
• 4. Select computing model (4) or (5) by user
• 5. Start calculation and visualization of each time step by calculating (4) or (5) and calculate the further computing grid for further

propagation step

• 6. User can export all results to frame and animation

Figure 9 shows an example of moving boundary of computing grid (red rectangle) at each time step. The domain is bounded by longitude 80E-105E and latitude 25N-10S (751x1051 pixels). The CPU time for calculation and real-time visualization of this example is about 1.9 mins on Laptop 1.4 GHz Pentium M Intel Processor. Figure 9. Tsunami propagation with moving boundary (resolution 601x610 pixels).

4.4. Results of numerical simulation

The tsunami arrived at the southern part of Phuket Island at first. It took 90 minutes to propagate from the tsunami source (Figure 10). Water surface descended first and ascended after half an hour. The tsunami arrival ran to the north direction along the coast of Thailand. Witnesses confirm this feature of the tsunami arrival time. The tsunami arrival time at Khao Lak was an hour later than on Phuket Island. Figure 10. Tsunami Time Line Every 30 minutes.

t=0.49 mins t=24.91 mins t=49.58 mins t=74.25 mins t=98.91 mins t=205.97 mins

• 5. CONCLUSION

The characteristics of the devastating Indian Ocean Asian Tsunami in various countries have made the understanding of the extreme event, in order to prepare early warning tools and systems which may allow an effective warning mechanism for the population, which was not in place or did not work effectively in the case of event. This paper introduces the tsunami propagation simulator model which calculated trends to check the model could predict the Tsunami arrival times and if it could be appropriate for a Disaster Warning System. In general it was found that the model can predict quite well the timing. The SiTProS as tsunami propagation simulator model has therefore been developed by ourselves, using a fast algorithm for a quick estimate of the Tsunami front propagation. We have tried to use of this model applied to the 26th December Tsunami and tsunami propagation in any areas, it shows a very accurate prediction of the time of the main events. The SiTProS Model runs in about 2.78 mins

• n a 2.8 GHz Intel Processor until the tsunami hitting Phuket based on domain grid 2101x2101 pixels.

In fact it was simulated that, in the case of 26th December 2004 Tsunami, the move of the originating point of a few degrees northward can cause a difference of about 1h in the predictions, due to bathymetry of the areas that in that particular area. The SiTProS model could be employed by experiment some points around the epicenter point to see its effect.

• 6. ACKNOWLEDGEMENTS

We would like to express our sincere gratitude and deep appreciation to Professor Fumihiko Imamura, DCRC, Tohoku University, Professor Yoshinobu Tsuji, ERI, University of Tokyo and Dr.Hiroyuki Matsumoto from JAMSTEC and Assoc. Professor Dr. Ahmet Cevdet Yalciner, Middle East Technical University, Civil Enigneering Department, Ocean Engineering Research Center, Ankara Turkey for their guidance, invaluable advice, supervision and material. We also heartily thank NOAA's National Geophysical Data Center (NGDC) for etopo2 data.

• 7. REFERENCES

[1] K. Abe, “Revised Mt and run-up estimate for the Indian Ocean Tsunami”, e-mail to ITIC Tsunami Bulletin Board posted on January 26, 2005. [2] L. Mansinha, and D. E. Smylie, “The displacement fields of inclined faults”, Bulletin of Seismological Society of America, vol. 61,

• no. 5, 1971, pp. 1433–1440.

[3] USGS : Largest Earthquakes in the world, http://neic.usgs.gov/neis/eqlists/10maps_world.html, referred on June 1, 2005. [4] T. Utsu, “Seismology”, pp. 376, Kyoritsu Shuppan, 2001. [5] H. Matsutomi, “Field Survey and Numerical Simulation on the 2004 Off-Sumatra Earthquake and Tsunami in Thailand”, http://www.tsunami.civil.tohoku.ac.jp/sumatra2004/C4.pdf, referred on June 1, 2005 [6] A. Annunziato and C. Best, “The Tsunami Event Analyses and Models”, the report at Institute for the Protection and Security of the Citizen Joint Research Centre European Commission, Jan 2005.