Numerical Study of Tsunami Propagation in Mentawai Islands West Sumatra

Shallow water equations were analyzed numerically for simulation of tsunami propagation from the source area to coastal zone of Mentawai Islands along West Sumatra Island. A triangular mesh was set for spatial discretization and a system of partial differential equation is reduced to a system of ordinary differential equations. Initial displacement based on Okada model was applied. Our numerical techniques were demonstrated with a simulation of tsunami generated from an earthquake with epicenter near South Pagai Island, Mentawai Islands, Indonesia.


Introduction
On the West Sumatra off coast, there are Mentawai Islands segment, which is a part of the Sumatra subduction zone, which results from the under-thrusting of the Indo-Australian plate beneath the overriding Eurasian plate, stretches over 5600 km [1]. In the last decade, these Mentawai Islands segment produced many great earthquake with potential tsunami events, including the October 25 th , 2010 Mentawai Islands tsunami earthquake. The earthquake occurred at 21:42:22.0 with the epicenter 3.487ºS, 100.082ºE and 20 km depth (NOAA Significant Earthquake Data), and generated a tsunami within estimated 8 minutes off the earthquake. Seismological analyses indicate that the Mentawai earthquake was a "tsunami earthquake", which produced a much larger tsunami than expected from the seismic magnitude [2].
Ulutaş [3] applied two kinds of tsunami numerical models based on the shallow water equations to simulate the South Pagai Island, Sumatra earthquake tsunami which occurred on 25 th , October 2010. The numerical models were analyzed by using nonlinear and dispersive long wave tsunami models, and earthquake parameters which are available from U.S Geological website. Satake, et. al [2] computed 25 th , October 2010, Mentawai Islands coastal tsunami heights from the tsunami source model which was proposed from a field survey and waveform recording with linear equations. The slip distributions on the fault planes were estimated through the inversion of tsunami waveform, whereas the strikes and rakes were estimated from USGS W phase solution.
In this study the propagation of tsunami waves generated in Mentawai Islands segment were studied numerically. A triangular mesh was employed for discretization of the governing equations, a system of partial differential equations, to a system of ordinary differential equations. The earthquake fault plane parameters based on Satake, et. al [2] were used for computation with Okada theory [4,5]. Numerical results illustrate the tsunami waves propagation generated by the 2010 Mentawai Islands earthquake.

Materials and Methods
A system of partial differential equations consisting of momentum equations and a continuity equation were used, and reduced by using discretization based on a triangular mesh to a system of ordinary differential equations.

Governing Equations
In the theory of long waves, the vertical acceleration of water particles are negligible compared to the gravitational acceleration except for an oceanic propagation of tsunami. Consequently, the vertical motion of water particles has no effect on the pressure distribution [6] which leads to two dimensional equations with dynamic and kinetic conditions, this is called the shallow water theory, can be obtained.
The numerical model is based on the nonlinear shallow water equations [7]. The following system of partial differential equations (1)-(3) were solve numerically for simulation of the tsunami propagation Variables and are the flux of -and -components, respectively which are defined by where and represent the -component and -component of velocity respectively, and = + ℎ is the total depth with sea depth ℎ and the water surface elevation from the mean sea level . The constant is the gravitational acceleration and is the Manning roughness coefficient.
System of equations (10)-(12) become the system equation (15)-(17) where = 2 ( +ℎ ) 7 3 The partial derivatives on the right-hand sides of the system equation (16)-(17) become The partial derivatives Values of partial derivatives at a nodal point are approximated by weighted averages of partial derivatives over elements that have the nodal point as vertices. For example, where 1 , 2 , … , are the areas of the elements which have the -th node as a common vertex, and � � ( ) is an approximate value of the partial derivative � � in the -th elements [7]. The collocation method was implemented in analysis of the system of equations above for discretization on a triangular mesh, and the system of partial differential equations was transformed to a system of ordinary differential equations. The system of those ordinary differential equations were solved numerically using a standard ODE solver i.e the fourth order Adam-Bashforth-Moulton predictor-corrector in conjunction with the fourth order Runge-Kutta method to generate values of approximate solutions at the first three steps.

Okada Model [4, 5]
Numerous theoretical formulations describing the deformation of an isotropic homogeneous semi-infinite medium have been developed with increasing completeness and generality of source type and geometry, and range from the derivation of the surface displacement due to a point source to the strain fields at depth [4]. The closed analytical expressions was checked and reviewed by Okada [4], which are already published to describe the surface deformation due to shear fault in a half-space. Furthermore, an unknown solution for the displacements strains and tilts arising from opening-mode dislocations.
Okada [5] established a model for calculation of internal deformation fields due to multiple sources that can be arbitrarily composed of shear and tensile faults of both point and rectangular types, where the surface deformation was formulated. This work is the extension of the previous work [4] to the internal deformation fields due to shear and tensile faults.
The Okada model was developed under the assumption that a rectangular fault plane is buried in a elastic semi-infinite half-space. On the case of sea surface deformation, the sea surface mimics the deformation of seafloor, since an earthquake usually occurs in seconds and the water column over the seafloor cannot escape within the short duration. The data of fault parameters which are required to compute the surface deformation include longitude and latitude, focal depth, length and width of fault planes, dislocation, strike angle, slip angle and dip angle.
Satake, et.al [2] divided the tsunami source area into 28 subfaults. The faults plane parameters [2] were used to computed the seafloor vertical deformation for a unit slip on each subfault.

Discussion
The system of partial differential equations (1)-(3) was analyzed for simulation of the tsunami propagation for Mentawai Islands case on October 25 th , 2010. Finite element analysis was applied to solve the system equations (1)-(3) numerically.

Bathymetry Data and Tsunami Source
The bathymetry data including latitude, longitude and depth from BDOC, National Oceanographic Database, GEBCO one minute grid [8] were used. Those data ranges from 85.0ºE to 106.0ºE and 15.0ºS to 6.0ºN. The data which originally given in form of latitude and longitude are transformed into projected coordinates by the Gauss-Krüger projection.
The fault plane is stretched from 99.86567ºE to 100.10976ºE and 2.55268ºS to 4.34144ºS which divided into 28 subfaults. Every subfaults have 30 km in length and 30 km in width. The strike and rake angle were 326º and 101º respectively, were estimated from USGS W phase solution. The dip angles were assumed to be 7.5º for 2-5.92 km depths subfaults and 12º for 9.83-16.07 km depths subfaults [2]. All those parameter values were employed in the Okada formula.
The vertical seafloor deformation for each subfaults which were computed using Okada model, was set as an initial condition to generate the initial tsunami wave. Figure  1 shows the initial surface deformation for the 28 subfaults in the fault plane source area. The results shown are the computational outcomes obtained by using Okada model. It shows the significant vertical surface displacement at several subfaults area, where the highest vertical deformation is approximately 3 m.

Simulation of the Tsunami Waves Propagation
The tsunami wave propagation is simulated for 1800 seconds simulation, which is equal to 1 hour tsunamis propagation. The initial tsunami wave form at 0 second, in the simulation area, is shown in Figure 2. Numerical results for tsunami wave propagation after 300 seconds and 600 seconds are shown in Fig. 3A (left) and Fig. 3B (right). Numerical results for the tsunami propagation after 900 seconds and 1200 seconds, 1500 seconds and 1800 seconds are shown in Fig. 4A (left) and Fig. 4B (right), Fig. 5A (left) and Fig. 5B (right), respectively.
In simulation initial tsunami wave form was ruptured before the 300 seconds, which is equal to 10 minutes after initial tsunami wave form generated. Satake, et. al [2] assumed that the tsunami arrival time was at 21:50 WIB on 25 th , October 2010. It is approximately 8 minutes after the earthquake which occurred at 21:42:22.0 WIB. The highest initial wave was approximately 3 m, and in simulation the highest wave of approximated 4 m was observed at several points in Mentawai Islands area, whereas the other wave heights are ranged from 1 m to 4 m. The results of this study showed the smaller wave heights compared with those reported on the field survey.
Field survey results of tsunami affected areas in the Mentawai Islands by Koresawa [9] reported that the highest and the lowest tsunami wave height at South Pagai Island was 9.3 m and 2.5 m, respectively. North Pagai Island was affected by tsunami wave height 8.8 m and 2.9 m as the highest and the lowest. Satake, et. al [2] indicated that the measured tsunami heights were mostly between 4 m and 7m, where the largest tsunami heights were measured on the central and southern South Pagai Island and south-western North Pagai Island.

Conclusions
Results from numerical study of tsunami event at Mentawai Islands, October 25 th , 2010, was presented. The Okada model was used to calculate the initial surface displacement, and the fault plane source parameters of the earthquake were utilized. A triangular element mesh was employed to discretize the governing equation.
The computational results of the initial surface deformation by Okada model shows that the highest initial wave height was 3 m. In simulation, tsunami wave reached the coastal areas of Pagai Islands in less than 300 seconds, at South Pagai Island and North Pagai Island, after the initial tsunami waves generated. In simulation the highest waves of approximately 4 m appear at several points around South Pagai Island and North Pagai Island. The simulation of this study showed the wave heights were smaller than those reported on the field survey. The tsunami waves were distributed along South Pagai Island, North Pagai Island, Sipura Island, edge of Siberut Island and reached in part of the coastal area of Sumatra Island at the province of Bengkulu.