SIMULATIVE ANALYSIS OF THE PROPAGATION OF A FAR-SOURCE HISTORICAL TSUNAMI ONTO PHILIPPINE COASTS

A simulative analysis methodology is presented and discussed to hindcast the propagation and shallow water transformation of a historical tsunami wave. The initial pulse of water surface induced is numerically modelled based on known geophysical data of earthquake magnitude and seismically induced seabed displacements. The propagation model accounts for the trans-sea movement, long wave propagation and damping, and shallow water transformations but excluding wave runup on the foreshore. The methodology is applied to the Philippine Trench 2012 tsunamigenic event using secondary data from regional geophysical databases and yielded good agreement with observed tsunami heights and arrival times recorded for local and regional locations, particularly at deeper and farther locations from the source.


INTRODUCTION
Tsunamis are potentially catastrophic natural hazards that are not widely studied in the Philippines. Considering the other calamitous hazards that affect the country in terms of economic loss, human casualty, damage to property and way of life, tsunamis are potentially catastrophic in terms of hazard intensity, spatial extent and resulting economic costs of damage. Research on tsunami propagation is relatively more extensive in the region, for example, in Japan, where highly catastrophic tsunamis have occurred more frequently in the past, leading to the early development of several studies on the science and engineering mitigating schemes for tsunamis. This includes tsunami forecasting as early as 1941, due to the occurrences of local tsunamis as early as the 1896 Meiji Great Sanriku Tsunami which claimed 22,000 lives with a run-up height of 38 m and the 1933 Showa Great Sanriku Tsunami which led to the development of tsunami mitigation schemes in the country (Shuto et.al, 2009). Up until today, strong tsunamis are frequently occurring in Japan, which is mainly due to their geographic location and crustal alignment. In 2011, a magnitude 9.0 was generated off the Pacific coast of Tohoku, Japan killing around 15,000 people and caused severe damage to a nuclear power plant, posing health risk concerns in the inundated area. In California, USA, several studies have also been conducted, i.e. Borrero et.al (2004) modelled locally generated tsunamis due to faulting and slope failures offshore of southern California using Okada's (1985) model for the initial condition and used the non-linear, shallow water wave equations to simulate propagation of the tsunami wave. The United States also conducts tsunami monitoring along their coastline and surrounding Pacific Coasts, which is managed by the National Oceanic and Atmospheric Administration -National Weather Service (NOAA-NWS) including studies on tsunami models on regional sources that may affect their coastline. In addition to this, one of the most catastrophic tsunami events is the Chile earthquake-generated tsunami of 1960, which showcased the trans-oceanic propagation of tsunami covering a wet distance of about 17,000 kilometers of ocean that caused movement of displaced seawaters as far as Izu, Japan on the other side of the Pacific Ocean. The Chilean tsunami of 1960 also led to the revision of design guidelines and upscaling of the design conditions for this well-known tsunami-prone area of the world. The South Asian Tsunami of December 2004 demonstrated the catastrophic magnitude and consequences of an infrequent tsunami event, which led to the loss of 300,000 human lives and damage to property of 1 billion US dollars.
In the Philippines, tsunamis are less regarded as the design basis of many engineering interventions for coastal protection and disaster mitigation, mainly due to the perceived low frequency of occurrence relative to the hazards of typhoons, extreme storm-induced waves, and coastal flooding. However, some research may also be found. One is a simplified model of wave incursion that is suited only for straight coastlines without blocking islands or complicated bathymetry (Cruz, 2010). A second related paper is on tsunami modelling in the West Philippine Sea, Dao et al. (2009), who conducted research on a hypothetical worst-case scenario tsunami in the West Philippine Sea due to movement of the Manila Trench. A generation and propagation model were carried out which showed results of wave heights of around 0.59 to 2.4 m upon reaching the Philippine coast. As there is a clear lack of tsunami studies and literature here in the Philippines, this study thus aims to attend to the following goals: (1) to develop a simulative methodology to study the generation, propagations and shallow water transformation of a historical tsunami; (2) incorporate the generation model into a coastal hydrodynamics model; and (3) to evaluate the applicability of the methodology by numerically hindcasting a far-source historical tsunami with observation data and is known to have reached Philippine coasts. The methodology is not intended to capture the tsunami wave runup and inundation processes on the foreshore slope and the hinterlands, which would require higher-resolution spatial and vertical data of both nearshore bathymetry and coastal topography.

PHILIPPINE TSUNAMI HAZARD
The Philippines is highly susceptible to tsunamis as the country is surrounded by active faults and trenches. The chances of tsunami generation from submarine and semi-terrestrial earthquakes are indeed high. Figure 1 shows the active faults and trenches in the Philippines. Two of these have historically produced significant tsunamis, namely, the Manila Trench which is located west of the island of Luzon, and the southern part of the Philippine Trench which bounds the islands of Visayas and Mindanao. Thus, it is important to determine tsunami susceptible areas to obtain the critical locations where more detailed analyses using prediction models is necessary. For the engineering interests, this capability is necessary in order to design the necessary mitigating measures prior to the occurrence of a catastrophe.

TSUNAMI SUSCEPTIBILITY
Tsunami susceptibility is related to the risk of a coastline to incur damage from the catastrophic effects generated by a tsunami. This damage is caused due to lack of protection of the coast as heightened tides and waves approach the coast. Since the country is surrounded by several active trenches which may cause tsunamigenic earthquakes at any time, our whole coastline is susceptible to the hazards generated by tsunamis. Thus, proper analysis prior to planning of mitigating schemes in case of a tsunami occurrence must be performed. Philippine tsunami susceptibility maps are currently present which give us a qualitative summary on the high-risk locations for a tsunami event

TSUNAMIGENIC EVENTS
Tsunami may be generated by any force that can produce a significant displacement of water that would create a series of large waves which propagate to the shore where damage to lives and property may occur. Tsunamis are usually generated from abrupt crustal movement from earthquakes, underwater landslides, and volcanic eruptions. In very rare cases, it may be caused by underwater blasts or explosions. As the Philippines is surrounded by active trenches, earthquakes events are the usual causes of tsunamis in the country. Several historical tsunamis have frequented the Philippine coastline. Figure 2 shows nearby recorded historical tsunami events along the southern coastline of Mindanao island.

TSUNAMI GENERATION BY SEISMIC SEABED DISPLACEMENT
Initial conditions of the model such as water surface displacement are governed by the initial seabed displacement at the tsunami source. The initial seabed displacement was modelled using the Okada (1985) model, which solves the seabed deformation due to shear and tensile faults in an elastic half-space assuming a homogenous and isotropic seabed. Given rectangular fault geometry (length, width, depth, strike, dip) and 3-component dislocation amplitude, shown in Figure 3, it computes the displacements, tilts, and strains at the free-surface. Equations for finite rectangular fault were used wherein fault dimensions were determined from the earthquake magnitude (Wells and Coppersmith, 1994). Then, the initial surface displacement was assumed equal to the displacement of the seabed and modelled as the initial water surface prior to simulation. The Philippine trench tsunami of 2012 was significant as it has reached in-bay tide gauges in Davao and Legaspi, which can supply the additional data comparison for model verification. However, runup heights at these locations are not high enough to be significantly remembered by residents. This lack of actual observation data from the localities of concern, prompted us to resort to water level measurements from gauges.

COASTAL ENGINEERING PROCEEDINGS 2020 4 GEOPHYSICAL DATA
Fault parameters for the 2012 Philippine earthquake event were obtained from the Global CMT Catalog (Ekström et.al, 2012). Which uses waveform data to determine earthquake source parameters (Dziewonski et.al, 1981). Two different characteristics of the fault movement are described in Table 1, both of which are simulated and compared with actual data. Fault length and width were estimated using a methodology by Papazachos et al. (2004) for a subduction zone dip-slip fault.

SECONDARY OBSERVATION DATA
Secondary observation data were obtained from the NOAA-NGDC, wherein observed tsunami wave heights at several locations were used as comparison to serve as verification of the propagation model for the tsunami event. These observation points vary from tide gauges for shallower areas and deep ocean gauges (Table 2).

SIMULATIVE ANALYSIS OF HISTORICAL TSUNAMI
Preliminary modelling data which consists of the deep crustal movement from the source, which then translates to the initial surface displacement, are necessary input prior to modelling of the tsunami propagation. This section discusses the model used to simulate the propagation of the 2012 Philippine trench tsunami.

COASTAL HYDRODYNAMICS MODEL
Tsunami propagation from a far-source earthquake was modelled using the hydrodynamics model MIKE 21, which is governed by the 2-dimensional depth averaged shallow water equations derived from the equations of conservation of mass and momentum. Bed roughness was considered with varying Manning's M ranging from 15 to 32 from shallower to deeper regions, respectively. Eddy viscosity was also considered, wherein horizontal eddy viscosity was modelled using Smagorinsky's formulation (1963). Consideration for Coriolis forcing is also included in the model. Computational domain and bathymetric modelling are further discussed in the succeeding sections.

Figure 5. Regional domain
The computational domain covers the entire area of the Philippines, extending out to the Philippine trench and further east to the corresponding deep ocean gauges in Guam and NE Manila. A finite difference spatial discretization scheme was used. Nesting was imposed on the shallower areas of Davao and Legaspi which are also enclosed by bays. The regional domain is shown in Figure 5. Grid nesting was performed with three nested grids to properly model dispersion and dissipation further into the bays up to the necessary tidal gauge locations. The smallest grid spacing used in the model was 90 m. The model discretization uses the Alternating Direction Implicit (ADI) technique and a double sweep algorithm. This leads to second-to third-order accuracy for the convective momentum terms and a well-conditioned solution algorithm (DHI, 2019). This is explained in more detail in the MIKE 21 Flow Model Scientific Documentation (DHI, 2019).
Ocean bathymetry was modelled using unstructured mesh with triangular finite elements, with gridded bathymetric data from GEBCO (2016) at 30 arc seconds, approximately 1 kilometer in spacing. Finer bathymetric data was additionally obtained for nested grids in shallower regions, these data were gathered from NAMRIA nautical charts at the areas of Legaspi and Davao. The finer resolution bathymetry for Davao Area is shown in Figure 6.

INITIAL OCEAN SURFACE DISPLACEMENT AND TSUNAMI TRACKING
The initial condition of the water surface was modelled using the Okada model (1985) wherein fault displacement was modelled with the fault parameters at rupture, i.e. depth, length, width, dip, slip, and strike. A subduction dip-slip model was assumed to obtain the length and width parameters from the earthquake magnitude (Papazachos et al., 2004). The Lame` coefficients which are the Shear modulus and lambda used were based on a four-layer half space model (Wang et al., 2003), which assumes these parameters to be dependent on the depth of rupture. The values used are 69.8 and 66.4 in units of giga-Pascals for the values of the Shear modulus and lambda, respectively. Both cases describing the fault movement were simulated, with parameters as shown in Table 1. A snapshot of the initial water surface is shown in Figure 7.

RESULTS AND DISCUSSION
Tsunami wave propagation is well observed from the model simulation results. Snapshots of the simulation succeeding the initial movement of the seabed are shown in Figure 8. Two surface waves COASTAL ENGINEERING PROCEEDINGS 2020 6 are observed to be generated in both east and west directions emanating from the tsunami source. The first wave of the tsunami quickly reached the eastern coastline of the Philippines in around 4 minutes. Snapshots of the propagation of the tsunami waves to the observation sites are shown for Legaspi, Davao, and the three deep water locations farther off the Pacific (Figure 8). A more detailed snapshot of the propagation of the tsunami wave to the observation sites housed inside the gulfs and bays in Davao and Legaspi (Figure 9). Higher water levels are observed to propagate into Legaspi Bay as compared to the Davao area, this is mainly due to the opening of the Legaspi coast to the east which results in minimal refraction as the wave approaches. However, both locations experience a reduced amplitude compared to coastlines directly under exposure, which is mainly driven by wave diffraction along the obstructing coastlines and islands. In the case of Davao, it was observed that tsunami crests were caused to bifurcate by Samal Island, located east of the observation point, and that the two separated waves eventually met inside the narrower channel, causing a heightened surge in water level in Davao. The model domain is interpolated based on the consolidated coarse to fine bathymetric data, into an unstructured mesh for numerical simulations, which extend past the Spratly archipelago where secondary wind data is obtained (Figure 8). The hydrodynamic model is based on the nonlinear shallow water equations (1,2) coupled with a spectral wave model derived from the wave action conservation equations (3,4). Considerations on wind friction, viscous effects, Coriolis forcing, bottom friction and depth-induced breaking are included in the coupled numerical model. Coupling of a hydrodynamic model with a wave model was deemed necessary in order to observe current-wave interaction at the project nearshore, which is highly varying given the changing waves and winds in both short-term and long-term seasonality. A similar approach was used by Bruneau et al. (2011) who modeled nearshore currents using a finite structured grid.

SYNTHESIS OF HISTORICAL TSUNAMI SIMULATION OUTPUTS
Based on the model simulation results, it is possible to calculate the tsunami travel times for the historical magnitude 7.6 Philippine trench tsunami of 2012. It is observed from the plots that the first wave of the tsunami first arrived at Homonhon Island, which is the easternmost part of Eastern Samar, at approximately 4 minutes after generation at the source. Also, the tsunami reached Guiuan, Eastern Samar at approximately 9 minutes after generation, propagating further into the entrance of Leyte Gulf 22 minutes after generation at the source. Siargao Island also experienced the surge, although not as severe as Eastern Samar. The first wave of the tsunami is seen to arrive at around 8 minutes along the north-eastern coastline of Siargao Island, Surigao del Norte.
This model was able to predict tsunami arrival times for deeper locations where tsunami wave dispersion is not as necessarily needed for an accurate model. Although higher order implementation shall produce more accurate results at certain locations. It is also necessary, to calibrate the model with several other Philippine tsunamis, of higher magnitude. After the confirmation of a properly calibrated model, such model may then be used to predict future tsunamis for several purposes of design of coastal protection structures and determine, for example, the proper set-back distances from the coastline in preparation for the occurrence of a catastrophic tsunami. This approach should be considered, as several lives and casualties, including economic loss and damages may be saved in the process of mitigating tsunami disasters.
In addition, studies on tsunamis in the Philippines should further be pushed for several reasons. This is driven by the fact that the Philippines is surrounded by several active faults and trenches, more studies and models need to be developed to properly predict runup heights, inundation extents, and arrival times of strong tsunamis that are bound to occur in the future. These results serve as necessary input to coastal protection, coastal management/planning, hazard mapping, establishing no-build zones, and setting easement requirements for no-build zone, and disaster response i.e. tsunami warning and evacuation. For this purpose, the availability of high-resolution data is of utmost importance to accurately predict inundation locations necessary for the above applications. Without these data, a tsunami numerical model will be inaccurate in predicting water level surges during a disaster.

CONCLUSIONS AND RECOMMENDATIONS
This paper studies generation and propagation of the 2012 Philippine Trench Tsunami. Data was compared at five (5) observation sites which showed good comparison at certain locations of the domain. From the simulation results, it is seen that observation points at shallower regions were inadequately modelled by the second-order equations of motion, notwithstanding the use of grid nesting and application of a depth-dependent Manning's roughness parameter. The dissipation of the tsunami wave is inadequately resolved and requires the addition of higher order terms to the model. On the other hand, the model is deemed to be generally adequate for other observation stations like Legaspi, Guam, and NE Manila which experienced less dispersion as the tsunami wave propagated to these stations. In order to properly model the far-source tsunami more accurately for certain locations, it is recommended to include higher-order terms in the model's governing equations.