Numerical Study on Tsunamis Excited by 2006 Pingtung Earthquake Doublet

The tsunami propagation and coastal run-ups generated by the 2006 Pingtung earthquake doublet is studied numerically. The numerical model COMCOT, which adopts a modified Leap-Frog finite difference scheme solving nonlinear Shallow Water Equations (NSWE), is employed. The assessment of the tsunami amplitude and extended effects is discussed in this paper using a set of nested bathymetric grids with a decreasing cell size focused on southern Taiwan. The three nested grids cover from the large far-field South-China Sea to the small near-shore regions. Several possible fault plane mechanisms are examined. Based on the arrival time at different tidal gauges, the suggested seafloor displacement is presented. Adopting the relocated seafloor movement with the fault parameters issued by Global CMT, the model results are favorably compared with the field tidal gauge wave form data. The result shows that the tsunami amplitude is about 35 cm and about 20-min period with the first depression at SyunGuaggZuei. A worse case scenario, Mw = 8.0, is studied. The maximum water levels can reach 4 m for Mw = 8.0 along the shoreline at HouBiHu. We also observe that owning to the source geometries and to the trapping of the tsunami energy by edge wave effect, the south-east coast of Taiwan coast is one of the most dangerous zones for the South China Sea earthquakes.


INTRODUCTION
An effective tsunami warning system hinges on accurate prediction of coastal run-ups immediately after the occurrence of a tsunamigenic earthquake.One way to achieve this is to establish a database of pre-calculated results for scenario earthquakes.Although major tsunamis are infrequent, numerical model can be validated by relatively moderate tsunami, i.e., the 26 Dec. 2006 Pingtung earthquake doublet.Among tidal gauge stations operated by the Central Weather Bureau, the most conspicuous tsunami signals resulting from the Pingtung earthquake doublet have been the observation at HouBiHu (Fig. 1 for location), with ~30 cm maximum peak-to-trough amplitude and ~20 min period for the first wave.The minor tsunami observed provides a way to test the effectiveness of simulation results.In this study, we compare the observed waves with tsunamis modeled to gain experience on how to establish an effective tsunami warning system in the future.
Firstly, we illustrate the application of the dislocation theory to obtain the static vertical displacements of ocean floor induced by the Pingtung earthquake doublet.All four possible combinations are considered due to the ambiguity of fault plane and auxiliary plane for each of the doublet.Secondly, the propagation of tsunami waves for the initial displacements is modeled using the COMCOT (see next section) software package, where the governing equations and methodology are introduced.Results show that the contribution of the second shock to the tsunami waves is relatively minor.Although the simulated waves are similar to those observed, it arrives about six minutes earlier than observed, which can be corrected by putting earthquake source 30 km further away from land.Upon simulation of a scenario M w 8.0 earthquake, the tsunami wave height at HouBiHu could reach four meters.

METHOD AND DATA
The process of a tsunami includes four steps: generation, propagation, amplification, and inundation.A detailed numerical model shall be able to simulate all the processes.In this study, the tsunami generation will be modeled based on an elastic model proposed by Manshinha and Smylie (1971) and a scaling law by Geller (1976).As for the tsunami propagation, amplification, and inundation, several numerical models can be used, such as the Japanese TIME model (Imamura et al. 1996;Imamura 1998), NOAA Center for Tsunami Research's MOST model (Titov and Synolakis 1998), and the Cornell COMCOT model.All these models are well validated, and the Cornell COMCOT model which solves nonlinear shallow water equations on the nested grids with the moving boundary algorithm will be adopted in this study.All the details of the models will be addressed later.We shall now introduce the data source and grid information.

Data Source and Grid Information
To simulation the tsunami propagation, three layers of grid are adopted.The largest layer with resolution of 2 minutes, referred to as Layer 1 (Fig. 2), covers the region between Hong Kong and Okinawa from west to east, and between Luzon and Suzhou from south to north.The second layer with 1 minute grid resolution is referred to as Layer 2 covers all the most interested area around Taiwan.The third layer (Fig. 3) with 0.125 minutes fine grid resolution covers the immediate epicenter region south of Taiwan.We shall address that all three layers are dynamically coupled to each other in the numerical models.

Generation
To model the vertical sea floor displacements caused by the Pingtung earthquake doublet, we adopt centroid moment tensor solutions obtained from the Global CMT Web Page (http://www .globalcmt.org).A detailed description of the CMT algorithm can be found in Dziewonski et al. (1981).Assuming a rigidity of 5 ´10 11 dyn cm -2 , a constant stress drop, and a unilateral slip on a rectangular fault plane where the length is twice as long as the width, we can infer the fault length, width and slip from the seismic moment using the scaling law (Geller 1976).Together with centroid depth and fault geometry described by strike, dip, and rake, the static vertical displacement of the ocean floor is modeled using the elastic dislocation theory (Manshinha and Smylie 1971;Okada 1985) (Fig. 4).We remark that the input source parameters (longitude, latitude, depth) represent the top center of the fault in the modeling scheme.We then follow the conventional approach to perform hydrodynamic simulation using the static deformation field as the initial condition.The approach is validated by the general rule that seismic rupture is much faster than water wave propagation.In practice, because there are two main shocks, each having two nodal planes as possible fault plane (Table 1), we use Case 12 to denote second nodal plane of the first shock.Case 11 and 12 are dominated by the negative displacement.The negative region has a dimension of 60 and 50 km for Case 11 and 12 (Fig. 5), respectively.The dimension is closely related to the wavelength of the generated tsunami.The largest negative displacement is -0.18 and -0.15 m for Case 11 and 12. Case 21 and 22 are thrust fault with relatively small amplitude of about 0.04 m for both Case 21 and 22 (Fig. 6).The combinations generate four sets of simulation results: Case 11+21; Case 12+21; Case 11+22; and Case 12+22.The simulated free-surface elevations from Case 11 and 12 at time = 474 sec will be superimposed with the surface elevations of Case 21 and 22.

Tsunami Propagation Model
The numerical method for the tsunami propagation and run-up, based on the shallow water approximation, has been adopted.This implies that the characteristic wavelengths of tsunami are considerably greater than the water depth and an averaging process over the water depth is performed.The appropriate mathematical model is the shallow-water theory.In the far-field and deep water region, the tsunami wave height is usually smaller then 1 m, which is much smaller than water depth.The advection terms can be neglected, and the linear shallow water wave equations read: where: x, y are the horizontal coordinates, h is the free-surface displacement, H = h + h is the total water depth, h is the still water depth, P = Hu, Q = Hv are the horizontal volume discharges, g is gravity, t is time, f is the Coriolis coefficient.When tsunami enters the near-shore region, the water depth reduces and the wave amplitude increases.The nonlinear advection and bottom friction terms become considerable and have to be included in the governing equations.The nonlinear shallow water wave equation reads: where t x and t y are the bottom frictions.
The bottom friction comes from Manning's formula and is expressed as: (2.5) where n is Manning's relative roughness.The Manning's relative roughness n is 0.013 in this paper to assume a sandy sea bottom.The Manning's n value has been tested in this study to have weak effect on the results.
In this study, the wave dispersion effect is assume to be small and can be neglected due to the near-field and long period tsunamis.However, for the cases of far-field tsunamis or wave of short period, the dispersion effect can play an important role.
The numerical method used in this study is COMCOT (Liu et al. 1994(Liu et al. , 1995;;Wang andLiu 2005, 2006) which is a tsunami modeling package.COMCOT adopts leap-frog finite differencing method to solve both linear and nonlinear shallow water equations in spherical and Cartesian coordinates.COMCOT has been widely used to investigate tsunami events, such as the 1992 Flores Islands (Indonesia) tsunami (Liu et al. 1994(Liu et al. , 1995)), the 2003 Algeria Tsunami (Wang and Liu 2005) and more recently the 2004 Indian Ocean tsunami (Wang and Liu 2006).With the flexible nested grids setup, we are able to perform the simulations from the near-coast region to the far-field region.In this study, linear shallow wave equations are solved in spherical coordinates on Layer 1 and 2. Nonlinear model is taken into account on Layer 3.
Boundary conditions have to be applied to solve for the governing equations.On the far-field open ocean, the radiation open boundary condition is adopted to permit the free wave propagation away from the numerical domain: (2.6) where v n is normal to the open boundary.The other type of the boundary adopted in this study is the internal boundary of the coastline, which located at the point on the boundary of the last sea points.In this paper, the moving boundary condition is adopted to consider the possible inundations and wave run-ups (Wang and Liu 2006).

RESULTS
Simulation results of tsunami generation and propagation are discussed in this section.The discussion will focus on the snapshots, comparison with tidal gauge data, time lag, epicenter location, maximum free-surface elevation profile, shoaling effect, and worst case scenarios.

Snapshots
Figure 7 shows the snapshots of free-surface elevation h on the Layer 2 of Case 11+22.The result shows that the first wave arriving in Taiwan is a negative wave and followed by a second large positive wave (times are 0 and 465 sec).The second earthquake happens at time = 474 sec.The effect on the second wave can be seen from the results between times 465 and 480 sec.The tsunami wave then enters the shallow area of Taiwan Strait with very low group velocity and relatively high wave amplitude, while the other part of the waves travels to the east coast of Taiwan due to the refraction (times 720 and 1080 sec).Because of the steep slope with deep water region, the wave height is small but the phase velocity is high (times 1800 and 7200 sec).In the shallow southwest coast of Taiwan, the high amplitude along the coastline from times 1080 to 7200 sec shows that some energy is trapped when the tsunami wave obliquely incident onshore.The trapped wave propagates along-shore at a much slower speed.Therefore, later at time 2550 sec wave crests in the deep sea are separated from the trapped crest.The along-shore propagation of the trapped tsunami waves can be easily identified by comparing the snapshots at time 2550 and 2850 sec.Note that the trapped tsunami lags behind the direct tsunami that propagates from the epicenter through the deep ocean to the coast.The existence of edge wave tsunami suggests that in a future tsunami incident the coastal zone should be prepared for incoming later tsunami arrivals even when the first few tsunamis have subsided.

Gauge Comparison
The tidal gauge data obtained from Central Weather Bureau (CWB) in Taiwan are compared with the numerical results.Data from four out of seven tidal gauge stations (Fig. 1) recorded significant tsunami signals will be examined in this paper.The recording sampling rate of the tidal gauge is 6 minutes while the computed digital series is every second.Figure 8 shows the gauge comparisons at tidal gauge stations JiangJyun, SyunGuaggZuei, HouBiHu, and Cheng-Gong.The tidal gauge data at SyunGuagGzue show that the crest of the tsunami arrives at time 19 minutes after the first earthquake occurrence.The largest wave amplitude is recorded at HouBiHu with 31.6 cm wave height at time 25 mi-  nutes.Because the tsunami arrival from the second earthquake is not obvious on the plot, the discussion shall be focused on the tsunami induced by the first earthquake.Inspecting the wave form, Case 12+21 and Case12+22 give closer wave form correlation with the tidal gauge data than Case 11+21 and 11+22.

Time Lag
Overall, the simulated results have similar wave period and wave height to those recorded by the tidal gauges.However, the time lag is about 6 minutes between the field measurements and the numerical results.In terms of the very near-field tsunami source, the 6-minute time lag is significant.This situation happens to all 4 simulations.It is difficult to attribute the problems to the entire tidal gauge stations, the pos-sible reasons are either the mis-location of the earthquake source or error in the assumed hydrodynamic model.Because the hydrodynamic model, COMCOT, is well calibrated in terms of the tsunami arrival time and the wave propagation, the distance between the closest gauge (SyunGuaggZuei) is short (less than 35 km), and dispersion has limited effect on the arrival time, we must assume that the reason is caused by the inaccuracy of the fault model and source location.Based on the phase velocity and the arrival time at the tidal gauges, data about the earthquake source should be re-estimated.The results from the re-located new epicenters will be discussed in the next section.

Relocate Epicenters
To be consistent with the tsunami arrivals, the epicenters of the Pingtung earthquake doublet have to be re-evaluated.
The epicenter of the first source is relocated at (21.36°N, 120.31°E), which is about 30 km to the south west of the epicenter estimated by Global CMT catalog.The new location of the second epicenter is at (21.76°N, 120.25°E).However, the correction on the second epicenter has less reliable owning to the unclear signal arrival shown on the time-history surface elevation.One piece of evidence to support the new surface displacement is the locations of the broken submarine communication cables.In this earthquake event, the submarine communication cables are broken at many locations (Fig. 9) (Data source: Chunghwa Telecom Fang-Shan Station).There are two reasons causing cable breaks.One is submarine mud flow, the other is the seafloor movement.
From the map of Fig. 9 we can see that the locations of the submarine communication cable breaks are around the relocated epicenter where the gradient of the surface displacement is large.Figure 10 shows comparisons of tidal gauge data and the numerical simulations with corrected epicenters.The comparison shows that the numerical results give better correlation with the observed data after relocating the epicenters.The results from Case 12+21 and 12+22 are closer to the tidal gauge data, and therefore Case 12+22 will be chosen as the one to be discussed.

The Maximum Free-Surface Elevation
Figure 11 shows the maximum free-surface elevation at each point from times 0 to 7140 sec. in Case 12+22.It reveals the information on the tsunami energy distribution and magnitude.The wave energy mainly travels towards the southern and southwest coast of Taiwan, where the water is shallow and the bathymetry is flat.The wave also travels along the steep east coast of Taiwan and attacks the northeast coast.

The Shoaling Effect
In the shallow water theory, the phase velocity C is estimated as C g h = + ( ) h .The dispersion relationship reads s 2 = gk tanh kh, in which s = 2p / T is the angular frequency, T is the wave period, k = 2p / l is the wave number, and l is the wavelength.When approaching the shoal, the water depth h decreases, wave velocity C slows down.However, the period T or angular frequency s remains approximately constant.Based on the dispersion relationship, if the wave number k increases, the wavelength l decreases.From the conservation of mass, the wave height a increases and the nonlinearity play a role.This shoaling effect, or the amplification effect, is significant while the long waves approach the shore.
Figure 11 shows the shoaling effect of the southern tip of Taiwan, where tsunami wave height is larger than other areas.During the shoaling process, the wave length reduces to about 15 km.It is roughly equal to the length of the flat shelf there.The shelf with a mild slope about 1/100 to 1/200 provides a great stage for tsunami to grow.The initial largest wave height is -0.18 m and followed by a similar bounced positive wave.The largest wave height reaches about 0.36 m along the shore.The amplification factor is about 2.0 at the off the southern tip of Taiwan.
Figure 12 shows the snapshot at time 900 sec and the cross-section profiles of sea surface elevation and bathymetry on the west side offshore of southern Taiwan.The shoaling effect takes place at the shelf area where the shelf length is close to the wave length.

The Worst Case Scenarios
In order to explore the potential tsunami hazard generated by large earthquake, a M w = 8.0 event was simulated.The fault parameters with the adjusted two epicenters are adopted to the elastic model.Based on the scaling law, the rupture length and width are 120.19 and 60.10 km, respectively.The dislocation is 3.83 m.The maximum wave height is recorded at HouBiHu (Fig. 13) at 4.0 m.The maximum wave height at SyunGuaggZuei is 2.2 m.We shall notice that at HouBiHu, the highest wave height does not occur at the first peak but at the second wave.The peak of the first wave arrives at HouBiHu at time 25 minutes, while the largest wave occurs at time 45 minutes.The wave period is about 20 minutes.
Figure 14 shows the maximum wave height of the case M w = 8.0.The dark red indicates the wave height is equal to or higher than 2.5 m.The map indicates that the most serious damage might happen to the southwest coast of Taiwan, where there are big cities with large population.The map also shows that 1.5 m wave height occurs at Taitung City and 1.0 m at Luodong Township.

CONCLUSIONS
The tsunami caused by the 26 Dec. 2006 Pingtung earthquake doublet is studied numerically.The numerical code, COMCOT, is adopted to perform the simulations.Three layers are coupled dynamically with one another to cover     The simulation results also show that the offshore regions of southern and southwestern Taiwan are vulnerable to tsunami attacks due to the 20 -30 km of gentle slope shelf bathymetry.Because of the steep slope and deep water, the tsunami energy is able to travel from the south to the north along the east coast with small energy loss to the northeast coast of Taiwan.
This paper also conducts a study on the worst case scenario.A M w = 8.0 event is simulated with the same fault parameter as that in Case 12 and modified epicenters.The result shows a maximum wave height of about 4.0 m at HouBiHu located at the southern tip of Taiwan.A map of maximum wave height for this worst case scenario is presented.

Fig. 1 .
Fig. 1.The bathymetry and tidal gauges of Layer 2. The red dots are two epicenters reported from Global CMT catalog.

Fig. 2 .
Fig. 2. The study area and the nested girds.The resolution of Layer 1 is 2 minutes; Layer 2 is 1 minutes; Layer 3 is 0.125 minutes.Linear shallow water equations are applied to Layer 1 and 2. Nonlinear shallow water equations are applied to Layer 3. Color indicates bathymetry.Black lines are the coastal lines.The red dots are the epicenters reported by Global CMT catalog.

Fig. 6 .
Fig. 6.The profiles of seafloor displacement/free-surface elevation of the second earthquake (Case 21 and 22).The color indicates the displacement.

Fig. 7 .
Fig. 7. (a) The snapshots of Case 11+22.The color indicates the free-surface elevation h.(b) The snapshots of Case 11+22.The color indicates the free-surface elevation h.(c) The snapshots of Case 11+22.The color indicates the free-surface elevation h.(d) The snapshots of Case 11+22.The color indicates the free-surface elevation h.(e) The snapshots of Case 11+22.The color indicates the free-surface elevation h.

2006
Fig. 8. (a) The time-history wave gauge measurement (blue lines) and numerical results.The epicenters are reported by Global CMT catalog.(b) The time-history wave gauge measurement (blue lines) and numerical results.The epicenters are reported by Global CMT catalog.

Fig. 9 .
Fig. 9.The suggested initial free-surface profile (or the seafloor movement).The red stars indicate the locations where submarine communication cable breaks (Data source: Chunghwa Telecom Fang-Shan Station).

Fig. 10 .
Fig. 10.(a) The time-history wave gauge measurement (blue lines) and numerical results of the re-located epicenters.(b) The time-history wave gauge measurement (blue lines) and numerical results of the re-located epicenters.
the off-shore to the near-shore of Taiwan.Both linear and nonlinear shallow water wave equations are solved on spherical and Cartesian coordinate systems for the linear and nonlinear models respectively.The effects of Pingtung earthquake doublet is simulated based on the wave height superposition.Four sets of faults parameters obtained from Global CMT catalog are simulated as Case 11+21, 11+22, 12+21, and 12+22.The simulated results are compared with the tidal gauge data from CWB.The comparison shows that Case 12 gives a better fit than Case 11.The sug-gested fault parameters for the first earthquake are: (Strike, Dip, Slip) = (329, 61, -98) and the focal depth = 19.6 km.There are no suggested fault parameters for the second earthquake.The simulation results also show that the arrival time of the wave crest is about 6 minutes faster than the tidal gauge measurements.This leads us to relocate the epicenter of the first earthquake to (21.36°N, 120.31°E).The results from the new epicenter show a better fit to the tidal gauge data in terms of the arrival time, maximum wave height, and the overall wave form.

Fig. 12 .
Fig. 12. (a)The distribution of free-surface elevation (upper plot), and the cross-section plot (lower plot).In the cross-section plot, both free-surface elevation (blue line) and the bathymetry (red line) are presented.(b) The distribution of free-surface elevation (upper plot), and the cross-section plot (lower plot).In the cross-section plot, both free-surface elevation (blue line) and the bathymetry (red line) are presented.

Fig. 13 .
Fig. 13.The time-history free-surface elevation from M w = 8.0 numerical simulation.Fig. 14.The maximum wave height from M w = 8.0 numerical simulation.