Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 1067002 10.3389/feart.2022.1067002 Earth Science Original Research An earthquake-triggered submarine mass failure mechanism for the 1994 Mindoro tsunami in the Philippines: Constraints from numerical modeling and submarine geomorphology Ramirez et al. 10.3389/feart.2022.1067002 Ramirez Alec Benjamin G. 1 * Ramos Noelynna T. 1 Nawanao Lyndon P. Jr 1 Mangahas-Flores Robelyn Z. 2 Narag Ishmael C. 2 Baba Toshitaka 3 Chikasada Naotaka 4 Satake Kenji 5 1 National Institute of Geological Sciences, College of Science, University of the Philippines Diliman, Quezon City, Philippines 2 Philippine Institute of Volcanology and Seismology, Quezon City, Philippines 3 Graduate School of Science and Technology, Tokushima University, Tokushima, Japan 4 National Research Institute for Earth Science and Disaster Resilience, Tsukuba, Japan 5 Earthquake Research Institute, The University of Tokyo, Tokyo, Japan

Edited by: Shaohong Xia, South China Sea Institute of Oceanology (CAS), China

Reviewed by: Polina Lemenkova, Université libre de Bruxelles, Belgium

Jean Roger, GNS Science, New Zealand

 *Correspondence: Alec Benjamin G. Ramirez, agramirez3@up.edu.ph

This article was submitted to Marine Geoscience, a section of the journal Frontiers in Earth Science

 30 11 2022 2022 10 1067002 11 10 2022 10 11 2022 Copyright © 2022 Ramirez, Ramos, Nawanao, Mangahas-Flores, Narag, Baba, Chikasada and Satake. 2022 Ramirez, Ramos, Nawanao, Mangahas-Flores, Narag, Baba, Chikasada and Satake

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Tsunamis have been known to result from a wide range of phenomena, such as earthquakes, volcanic eruptions, submarine mass failures, and meteorite impacts. Of earthquake-generated tsunamis, those arising from strike-slip mechanisms are less common, with the 1994 Mindoro tsunami in the Philippines among the few known examples. The 1994 Mindoro tsunami followed a Mw 7.1 earthquake along the right-lateral Aglubang River Fault. The tsunami affected the coasts surrounding the Verde Island Passage, one of the Philippines’ insular seas located between the islands of Luzon and Mindoro, and east of the West Philippine Sea margin. A total of 78 lives were lost due to the earthquake and tsunami, with 41 being directly attributed to the tsunami alone. Despite the close spatial and temporal association between the 1994 Mindoro earthquake and tsunami, previous numerical modeling suggests the need for other contributing mechanisms for the 1994 tsunami. In this study, we conducted submarine geomorphological mapping of the South Pass within the Verde Island Passage, with particular focus on identifying possible submarine mass failures. Identification of submarine features were based on Red Relief Image Map (RIMM), Topographic Position Index (topographic position index)-based landform classification, and profile and plan curvatures derived from high-resolution bathymetry data. Among the important submarine features mapped include the San Andres submarine mass failure (SASMF). The San Andres submarine mass failure has an estimated volume of 0.0483 km3 and is located within the Malaylay Submarine Canyon System in the Verde Island Passage, ∼1 km offshore of San Andres in Baco, Oriental Mindoro. We also explored two tsunami models (EQ-only and EQ+SMF) for the 1994 Mindoro tsunami using JAGURS. The source mechanisms for both models included an earthquake component based on the Mw 7.1 earthquake, while the EQ+SMF also included an additional submarine mass failure component based on the mapped San Andres submarine mass failure. Modeled wave heights from the EQ-only model drastically underestimates the observed wave heights for the 1994 Mindoro tsunami. In contrast, the EQ+SMF model tsunami wave height estimates were closer to the observed data. As such, we propose an earthquake-triggered, submarine mass failure source mechanism for the 1994 Mindoro tsunami.

 1994 Mindoro earthquake and tsunami submarine geomorphology submarine mass failure numerical modeling Verde Passage Philippines 1 Introduction

Tsunamis have been a ubiquitous phenomenon around the world, with several notable examples occurring within the past 2 decades. While most destructive tsunamis are generated by major to great megathrust earthquakes, such as the 2004 Indian Ocean and the 2011 Tohoku-oki tsunamis, more recent events arose from less common, non-seismic sources. These include the 2022 Hunga Tonga tsunami, which followed the eruption of the Hunga Tonga-Hunga Ha’apai submarine volcano, and the 2018 Sulawesi tsunami, which resulted from the combined effects of a Mw 7.5 earthquake along the Palu-Koro fault system in Indonesia (Heidarzadeh et al., 2018), associated submarine landslides (Gusman et al., 2019; Pakoksung et al., 2019), and coastal liquefaction (Sassa and Takagawa, 2018).

With exception to impact-generated tsunamis, tsunamis from strike-slip earthquakes are even much less common, as the resulting vertical seafloor displacements from such events are usually too small to initiate tsunamis (Tanioka and Satake, 1996; Heidarzadeh et al., 2017). Aside from the previously mentioned 2018 Sulawesi tsunami, the 1994 Mindoro tsunami in the Philippines is another example of a tsunami occurring after a strike-slip earthquake (Tanioka and Satake, 1996; Pakoksung et al., 2019).

The 1994 Mindoro tsunami was preceded by a Mw 7.1 earthquake along the Aglubang River Fault, within the South Pass of the Verde Island Passage in the Philippines. The tsunami produced wave heights at shore exceeding 7 m (Supplementary Table S1), despite the earthquake’s dominantly right-lateral, strike-slip mechanism. Tanioka and Satake (1996) suggested that the coseismic horizontal displacement of the coastal shelf, which is intersected obliquely by the right-lateral Aglubang River Fault, played an important role in producing the observed tsunami wave amplitudes. Contrary to this, numerical modeling by the same authors has shown that a source mechanism based on Mw 7.1 earthquake alone severely underestimates the wave height amplitudes and distribution compared to actual observations of the tsunami run-up and inundation.

The possible role of a submarine mass failure (SMF) in the generation of the 1994 Mindoro tsunami was then proposed by Tanioka and Satake (1996), and is explored in detail in this study. In general, SMF-generated tsunamis are characterized by much larger run-ups in the near-field compared to the far-field (Okal and Synolakis, 2004), which was also observed for the 1994 Mindoro tsunami (Imamura et al., 1995). Moreover, SMFs have been invoked as an important contributing mechanism for tsunamis similar to the 1994 Mindoro tsunami, wherein the estimated sea surface displacement derived from the associated earthquake is insufficient to account for the resulting run-ups. These include the 1946 Alaska (Fryer et al., 2004; von Huene et al., 2014), the 1998 Papua New Guinea (Synolakis et al., 2002; Satake and Tanioka, 2003), and the 2018 Sulawesi (Gusman et al., 2019; Pakoksung et al., 2019) tsunamis.

 2 Philippine tectonics and tsunamis

The archipelagic nature of the Philippines, combined with its tectonic setting, makes it susceptible to tsunamis (Bautista et al., 2012). Most of the Philippine archipelago is within the Philippine Mobile Belt (PMB), which is bounded by oppositely-dipping subduction zones to its east and west, with the exception of the Palawan Microcontinental Block (PMCB) (Figure 1A; Gervasio, 1971; Rangin, 1991; Yumul et al., 2008; 2009; Lagmay et al., 2009; Mines and Geosciences Bureau, 2010). To the east, the Philippine Sea Plate (PSP) is subducting beneath the PMB along the Philippine Trench (Hamburger et al., 1983; Ozawa et al., 2004; Yumul et al., 2008; Mines and Geosciences Bureau, 2010). To the west, the oceanic lithospheres underlying the South China Sea, Sulu Sea, and Celebes Sea are subducting beneath the PMB along the Manila, Negros-Sulu, and Celebes trenches, respectively. (Hayes and Lewis, 1984; Mitchell et al., 1986; Rangin et al., 1999; Yumul et al., 2008; Mines and Geosciences Bureau, 2010). Of the 35 certain tsunami events in the Philippines with wave heights of at least 1 m (Table 1; Bautista et al., 2012; NCEI/WDS, 2022), three of the largest were generated by subduction-zone earthquakes. These are the 1897 Zamboanga tsunami, which was initiated along the Sulu Trench, and the 1918 Celebes Sea and 1976 Moro Gulf tsunamis which were initiated along the Cotabato Trench. The 1976 Moro Gulf tsunami (Badillo and Astilla, 1978; Bautista et al., 2012; Claro et al., 2021) is the most destructive of the known tsunami events in the Philippines.

Philippine tectonic setting and notable tsunamis. (A) The 35 positive tsunamigenic earthquakes in the Philippines since 1589 AD with recorded or estimated wave heights at shore of at least 1 m (Bautista et al., 2012; NCEI/WDS, 2022) are plotted with yellow circles, or with the corresponding focal mechanism solution (FMS) diagrams. The most significant tsunamis, with maximum wave heights at shore exceeding 4 m, are also labeled. Oceanic trenches are also shown and labeled, including active faults in red (PHIVOLCS). (B) Location of the South Pass within the Verde Island Passage, where the 1994 Mindoro tsunami occurred. (C) Tsunami wave heights (modified from Imamura et al., 1995; see Supplementary Table S1) and directions of wave inundation (PHIVOLCS, 1994). The trace of the offshore extension of the Aglubang River Fault is also indicated, as mapped from the analysis of high-resolution bathymetry data in this study. FMS diagrams are from the Global Centroid Moment Tensor catalog (GCMT; Dziewonski et al., 1981; Ekstrom et al., 2012). PMB = Philippine Mobile Belt; PMCB = Palawan Microcontinental Block; PSP = Philippine Sea Plate; MT = Manila Trench; NT = Negros Trench; ST = Sulu Trench; CT = Cotabato Trench; ELT = East Luzon Trough; PT = Philippine Trench; PF = Philippine Fault; SS = Sulu Sea; CS = Celebes Sea; KIG = Kalayaan Island Group; PR = Philippine Rise; ARF = Aglubang River Fault; CMF = Central Mindoro Fault; LVPF = Lubang-Verde Passage Fault; BR = Baco River; MR = Malaylay River.

The 35 tsunamigenic earthquakes in the Philippines, with recorded or estimated maximum wave heights at shore of at least 1 m.

Id Y M D Hr Min Long Lat Depth (km) Mag a Wave height b Places affected
1 1828 11 9 18 30 119.5 13.7 - 6.6 1 Port of Manila
2 1840 3 22 8 30 123.85 12.95 - 6.8 2 Sorsogon Bay
3 1863 6 3 19 30 120.9 14.55 - 6.5 1–2 Manila Bay
4 1869 8 16 15 0 123.85 12.4 - 6.5 1 Masbate, Masbate
5 1880 7 18 12 40 121.55 14.9 - 7.6 2 Puerto Real, Quinanliman and Tacligan in Real, Quezon
6 1897 9 21 13 12 121.3 7.5 - 7.5 6 Sulu Sea
7 1917 1 31 12 2 125.6 5.5 - 6.4 1.5 Glan, Sarangani
8 1918 8 15 20 18 124.4 5.9 - 8 8 Celebes Sea coast; Lebak, Sultan Kudarat
9 1921 11 12 2 36 127 8 - 7.5 2 Manay, Mati, and Caraga in Davao Oriental
10 1923 7 18 10 42 125 9.3 - 5.5 1 Mambajao, Camiguin
11 1924 4 15 0 20 126.5 6.5 - 8.3 2 Pujada Bay
12 1924 5 6 0 16 119 16 - 7 1 Agno, Pangasinan
13 1925 5 5 18 7 122.7 9.3 - 6.8 1–2 Southern coast of Negros
14 1925 5 25 11 44 122.5 12.5 - 6.25 2 Tugdan, Romblon
15 1925 11 13 20 14 125 13 - 7.3 2 Batag Island, Northern Samar
16 1928 6 15 14 13 121.5 12.5 - 7 1 Mangarin Bay
17 1928 12 19 19 27 124 7 - 7.5 1 Cotabato River and Illana Bay
18 1929 6 13 17 24 127 8.5 - 7.2 2 Hinatuan Bay
19 1934 2 14 11 59 119 17.5 - 7.6 2 San Esteban and Vigan, Ilocos Sur
20 1937 8 20 19 59 121.5 14.5 - 7.5 1 Lopez and Calauag Bays
21 1948 1 25 1 46 122 10.5 - 8.2 2 San Joaquin, Miagao, Oton in Iloilo; Nueva Valencia, Guimaras
22 1968 8 2 4 19 122.2 16.5 - 7.3 3 Casiguran Bay
24 1970 4 7 13 34 121.717 15.761 - 7 4.3 Baler and San Luis, Aurora
25 1973 3 17 16 30 122.787 13.372 - 7.3 1.3 Calauag and Alabat, Quezon
26 1975 10 31 16 28 125.993 12.54 - 7.4 1–2 Taft, Eastern Samar
27 1976 8 17 0 11 124.023 6.262 33 8.1 9 Lebak, Sultan Kudarat; Moro Gulf
28 1990 2 8 15 15 124.694 9.755 16.2 6.6 1 Alijuan River, Duero, Bohol; Bohol Strait; Camiguin Island
29 1990 7 16 16 26 121.172 15.679 15 7.8 2 Brgy. Darigayos, Luna, La Union
30 1992 5 17 17 49 126.753 7.183 34 7.1 5 Manay and Caraga, Davao Oriental
31 1994 11 15 3 15 121.087 13.532 15 7.1 7 Northern Oriental Mindoro; Verde Island in Batangas
32 1995 4 21 8 34 125.58 12.059 21.7 7.3 2 San Julian and Dolores, Eastern Samar
33 1999 12 12 2 3 119.67 15.85 35.1 6.8 1.5 Iba and Palauig, Zambales
34 2002 3 6 5 15 124 6.1 28.7 7.2 2 Palimbang, Sultan Kudarat
35 2012 2 6 11 49 123.14 9.97 12 6.9 2 La Libertad, Negros Oriental

aMoment magnitude. For earthquake magnitudes reported as surface wave (Ms) and body wave (Mb) magnitudes, empirical equationsby Scordilis (2006) were used to convert the values to moment magnitude.

Maximum wave height at shore (meters).

The row in bold emphasizes the 1994 Mindoro tsunami, which is the focus of this study.

Earthquakes along offshore intraplate and upper-plate faults have also produced tsunamis in the Philippines—the most notable of which are the 1970 Baler, 1992 Davao Oriental, and 1994 Mindoro tsunamis (Bautista et al., 2012). The 1994 Mindoro tsunami, which is the focus of this study, occurred within the Verde Island Passage (Figure 1B). Major tectonic structures in proximity to the Verde Island Passage include the southern terminus of the Manila Trench, and active strike-slip faults such as the WNW-trending left-lateral Lubang-Verde Passage Fault (Karig, 1983; Rangin et al., 1988), the NNW-trending Central Mindoro Fault (Karig, 1983; Rangin et al., 1988), and the NNW-trending right-lateral Aglubang River Fault (PHIVOLCS, 1994).

 3 The 1994 Mindoro tsunami

In the early morning of 15 November 1994, at around 3:15 a.m. local time (GMT+8), a Mw 7.1 earthquake occurred underneath the South Pass of the Verde Island Passage in the Philippines (Figures 1B–D; PHIVOLCS, 1994). The earthquake is associated with the offshore extension of the right-lateral Aglubang River Fault, based on its epicenter (13.5 N, 121.1 E), focal mechanism solution, and surface rupture in Oriental Mindoro along the Aglubang River Fault (PHIVOLCS, 1994). Within a few minutes after the earthquake, tsunami waves hit the surrounding coastal areas of northern Oriental Mindoro, Baco Islands, Verde Island, and Batangas in Southern Luzon (PHIVOLCS, 1994; Imamura et al., 1995). Among the hardest hit is Barangay San Andres in the town of Baco, Oriental Mindoro, and the Baco Islands, which were inundated by tsunami waves exceeding 7 m and 5 m, respectively (Supplementary Table S1). A total of 78 lives were lost, with 41 directly due to the tsunami. The cost of repairing and rehabilitating the damages to infrastructure was estimated at around PHP 5.5 million (PHIVOLCS, 1994). The 1994 tsunami also left behind a 1-m high tsunami beach ridge along the coast of Barangay San Andres (PHIVOLCS, 1994), although it was not observed during fieldwork conducted in 2019 (Figure 2). The tsunami beach ridge was most likely eroded given that the Barangay San Andres coast has been subjected to at least 120 m of coastal erosion since 1992, as observed from multiple Landsat and Sentinel-2 satellite imagery (Supplementary Figure S1). An extreme wave event deposit, most likely a preserved section of the 1994 Mindoro tsunami, was also identified in a geoslicer sample obtained in Barangay San Andres (Figure 2C). The suspected 1994 tsunami deposit is characterized by the presence of large (>5 cm) mud clasts in a muddy, medium to coarse sand matrix, with a basal erosional contact with the underlying, mangrove swamp sediments.

Barangay San Andres in Oriental Mindoro, 25 years after the 1994 tsunami. (A) NNE-looking image overlooking Sitio Malaylay in Barangay San Andres, taken during low tide. The area is mostly a mangrove swamp, with Nypa fruticans as the most common mangrove. Note the dilapidated structures which were initially damaged by the 1994 tsunami. The location of Figure 2A is shown in Figure 1C. The exposed tidal bar in the foreground of B is also labeled. (B) SE-looking, along-the-coast view of Sitio Malaylay. Note the protruding trunks of drowned mangrove trees along the foreshore and the absence of the tsunami beach ridge that was deposited by the 1994 tsunami, both of which indicate post-1994 coastal erosion. The location of a geoslicer sample, MAL-G3, is also shown. (C) Visual stratigraphic interpretation of MAL-G3. The 1994 Mindoro tsunami deposit is partially preserved, overlying the older mangrove swamp sediments.

The 1994 Mindoro tsunami was initially assumed to be caused by the Mw 7.1 earthquake (PHIVOLCS, 1994; Bautista et al., 2012), due to the close spatial and temporal association between the two. However, as previously mentioned, results of previous numerical modeling of the 1994 tsunami based on the Mw 7.1 earthquake alone (Tanioka and Satake, 1996) are inconsistent with observed tsunami wave heights. As such, in this study, we explored the possibility of one or several submarine mass failures as a contributing mechanism for the 1994 tsunami.

 4 Submarine geomorphology of the South Pass 4.1 Morphometric analysis and geomorphological interpretation of bathymetry data

The 20 m resolution digital bathymetry data for the Verde Island Passage (Supplementary Figure S2), obtained from multibeam surveys conducted by the National Mapping Resource and Information Authority (NAMRIA) in 2011, was used in the submarine geomorphological mapping. Geomorphometric analysis of the bathymetry data was performed through topographic position index (TPI)-based landform classification and derivation of land surface parameters (i.e., slope, openness, profile and plan curvature) in a GIS platform. These parameters were then used as basis in identifying submarine geomorphological features (i.e., gullies, ridges, fault line scarps, shelf breaks and scarps, seafloor bedforms, and submarine mass failures).

4.1.1 Delineation of linear features using red relief image maps

The Red Relief Image Map (RIMM; Chiba et al., 2008) is a pseudo-3D relief visualization technique that highlights both coarse- and fine-resolution surface features simultaneously across a wide variety of topographic positions, without the dependency on light directions, unlike shaded relief maps. It is generated by combining two land surface parameters: slope and differential openness. In RIMM, a red color ramp is used to visualize the slope angle, while its saturation is modified based on the differential openness. As such, convex features (e.g., ridges and crests) appear as bright areas; concave features (e.g., valleys and depressions) and steep slopes appear as red; and flat areas appear as gray.

The RRIM for the South Pass was generated in QGIS 3.22.4. The slope was calculated using the Slope terrain analysis module (GDAL/OGR contributors, 2020), following Horn’s algorithm. The positive (PO) and negative openness (NO) were calculated using the Topographic Openness module from the System for Automated Geoscientific Analyses (SAGA) (Conrad et al., 2015), following the definition of openness by Yokoyama and others (2002). We then calculated the differential openness (DO) from the PO and NO, using the formula from Chiba et al. (2008): D O = ( P O − N O ) 2

The slope was visualized using a red color ramp, while the differential openness used a black-to-white color ramp. The RRIM effect is achieved by overlaying the slope over the differential openness, using the blending mode Multiply. Prominent linear features were manually delineated from the RRIM (Figure 3A).

DEM-derived layers used in the submarine geomorphological mapping of the South Pass, Verde Island Passage. (A) Red Relief Image Map (RRIM) with manually delineated linear features. (B) Landform classification based on Topographic Position Index. (C) Profile curvature. (D) Plan curvature.

 4.1.2 TPI-based landform classification

The topographic position index (TPI) is calculated by obtaining the difference between the elevation of the point of interest (Z0) and the mean elevation (ZC) around a neighborhood area with a predetermined size C (Guisan et al., 1999; Wilson and Gallant, 2000). By calculating the TPI using two neighborhood sizes (CS and CL), distinct small- and large-scale landform units can be defined based on the values of TPI and the standard deviation (SD) within C (Weiss, 2001; De Reu et al., 2013).

The TPI Based Landform Classification module in SAGA (Conrad et al., 2015) was used to derive the landform classes, using pixel (px), small (CS), and large neighborhood (CL) sizes of 20 m, 100 px, and 1000 px, respectively. The landform classes for submarine geomorphology were modified from those originally defined by Weiss (2001) for terrestrial landscapes (Table 2).

Landform classification based on TPI values.

Small neighborhood (CS = 100) Large neighborhood (CL = 1,000) Landform classes Id
Weiss (2001) This study
Z0 < -SD Z0 < -SD Canyons, deeply incised streams Gullies, local incised valleys and depressions 1
Z0 < -SD 0 ≤ Z0 ≤ SD Midslope drainages, shallow valleys
Z0 < -SD Z0 > SD Upland drainages, headwaters
-SD ≤ Z0 ≤ SD Z0 < -SD U-shaped valleys Wide valleys, canyons, troughs, and depressions 2
-SD ≤ Z0 ≤ SD 0 ≤ Z0 ≤ SD; slope ≤5 Plains Plains, plateaus, and shelves 3
-SD ≤ Z0 ≤ SD 0 ≤ Z0 ≤ SD; slope >5 Open slopes Open slopes 4
-SD ≤ Z0 ≤ SD Z0 > SD Upper slopes Upper slopes 5
Z0 > SD Z0 < SD Local ridges/hills in valleys Midslope and local ridges 6
Z0 > SD 0 ≤ Z0 ≤ SD Midslope ridges
Z0 > SD Z0 > SD Mountain tops, high ridges Seamounts, high ridges, and scarps 7

The modifications include the use of appropriate submarine landform terms corresponding to each landform class. For example, additional landforms such as shelves and plateaus were added to the submarine landform class 3, since both are similar to plains in terms of their low relief. Moreover, Weiss’ landform classes 1, 2, and 3 were merged into a single class grouping together all narrow channel-like landforms (i.e., gullies, upland drainages, and canyons), as well as narrow valleys and depressions. This is to simplify the distinction of negative relief landforms into narrow (class 1) and wide (class 2) morphologies. Similarly, Weiss’ landform classes 8 and 9 were also merged into a single class, which includes both midlsope ridges and local ridges located within valleys. This again, allows for a simpler distinction versus prominent ridges (high ridges, or class 7 in this study).

The TPI-based landform classification for the South Pass is shown in Figure 3B. Violet and blue were used for landform classes with negative relief (1 and 2); yellow for low relief landforms (3); orange for open and upper slopes (4 and 5); and green for positive relief landforms (6 and 7).

 4.1.3 Profile and plan curvatures

Profile and plan curvatures are measures of the rate of change (radians per m) of the slope and aspect of a topographic surface, respectively. Profile curvature is measured along the vertical direction. On the other hand, plan curvature is measured along the horizontal direction, parallel to elevation contours, and hence are also referred to as contour curvature. Positive and negative values for both plan and profile curvatures correspond to convex and concave morphologies in relation to the downslope direction (Wilson and Gallant, 2001). The profile and plan curvatures in 20-m resolution for the South Pass were also calculated in 20-m resolution using SAGA, as shown in Figures 3C,D.

 4.1.4 Geomorphological classification of linear features

The linear features delineated from the analysis of RRIM for the South Pass were manually classified based on their interpreted geomorphology. The TPI-based landform classification, profile curvature, and plan curvature were also used as guides for interpretation. A summary of the submarine geomorphological features and their expressions in the RRIM, TPI-based landform classification, and profile and plan curvatures are shown in Table 3. The resulting submarine geomorphological map and bathymetric profiles of the South Pass are also shown in Figure 4.

Geomorphological classification of linear features.

Landform Description a TPI landform class b Profile curvature Plan curvature
Gullies Narrow, incised channels; flanked by midslope ridges 1 concave
Ridges Linear bathymetric highs flanked by steep slopes 6, 7 convex convex
Fault line scarp Linear slope breaks or steep valleys, conforms with the focal mechanism of the 15 November 1994 Mw 7.1 Mindoro earthquake c 1 concave -
Shelf breaks, scarps Linear slope breaks separating an upper, gentler slope face or shelf and a lower, steeper slope face 7; separates 3, 5 from 4 convex -
Seafloor bedforms Low relief, regularly-spaced crests and troughs 3 - -
Submarine mass failures d
Headwall scarps Steep, arcuate and concave-facing downslope 7 convex concave
Lateral margins Scarps or ridges perpendicular to the headwall scarp 6, 7 convex
Toe deposits High surface roughness due to displaced blocks; occurrence of pressure ridges

As recognized from the red relief image map (RRIM).

See Table 2 for the TPI, landform categories.

Oriented subparallel to nodal plane one from GCMT, focal mechanism (Strike: N 21 W; Dip: 70; Rake: 178; Dziewonski et al., 1981; Ekstrom et al., 2012).

From morphological descriptions of submarine mass failure structures by Scarselli (2020) and morphological criteria for the identification of terrestrial landslides from high-resolution digital terrain data by Van Den Eeckhaut and others (2012) and Pawluszek (2019).

Submarine geomorphology of the South Pass, Verde Island Passage. (A) Map of submarine landforms and features. The headwall scarps of the 15 mapped probable SMFs (Table 4) are also labeled. The extent of the 3D rendering in Figure 5 is also shown. (B) Relief map of the South Pass, labeled with identified main submarine features and transects for the bathymetric profiles shown in (C). ARF = Aglubang River Fault; BIS = Baco Islands Shelf; CS = Calapan Saddle; MSCS = Malaylay Submarine Canyon System; NMS = Northern Mindoro Shelf; VBIS = Verde-Baco Islands Saddle; VIS = Verde Island Shelf; VPB = Verde Passage Basin; BR = Baco River; MR = Malaylay River.

 4.1.5 Mapping of submarine mass failures (SMFs)

Particular focus was also given on the identification of probable SMFs which were primarily based on their morphological expressions as discussed by Scarselli (2020), such as headwall scarps, lateral margins, and toe deposits. Morphological criteria for the identification of terrestrial landslides from high-resolution digital terrain data (Van Den Eeckhaut et al., 2012; Pawluszek, 2018) were also used to supplement the identification of SMFs, based on the morphological similarities between submarine and terrestrial landslides. The head scarp length, transport direction, and area of each mapped SMF was also measured in QGIS.

Furthermore, we also identified the SMFs that have probably occurred after the 1950s, based on the comparison between the 2011 bathymetry with 1950s point depth sounding data. The point depth sounding data was obtained from NAMRIA topographic map of Calapan, Oriental Mindoro (3260 III), which was based on aerial photography and coastal surveys and offshore surveys from 1947 to 1953.

Both bathymetry data sets generally show minor differences with each other, with a coefficient of determination (R 2) of 0.95 (Supplementary Figure S3). As such, we can use the two bathymetry datasets to reliably estimate the net erosion from 1956 to 2011. The net erosion for each point depth location was estimated based on the difference between the two datasets. Post-1956 SMFs were identified based on the occurrence of net erosion downslope of SMF scarps.

 4.2 Main undersea features

The submarine geomorphological map and bathymetric profiles of the South Pass is shown in Figure 4. The South Pass (SP) is divided into two regions of distinct bathymetric levels by a ∼21 km long SW-facing slope, trending NW from the Northern Mindoro Shelf (NMS) towards Verde Island. The Verde Passage Basin (VPB) floors the SW region and further extends outside SP to the northwest, with its shallowest portions at the 300-m isobath. To its south, the VBP is bounded by the western Northern Mindoro Shelf (wNMS). The wNMS is composed of a gentle (<5°) upper shelf with a maximum width of ∼3 km, and a steeper (5–10°) insular slope transitioning into the VPB.

On the other hand, the NE region of SP is much shallower than VPB, and is divided into a series of insular shelves and saddles. The eastern Northern Mindoro Shelf (eNMS) is also composed of a gentle (<5°) upper shelf underlying the Calapan Bay, with a maximum width of ∼1.3 km, and a steeper insular slope extending to the 150-m isobath, where it transitions to the Calapan Saddle (CS). CS is ∼9-km long, ∼1.2-km wide, and has a maximum depth of less than 250 m. To its north lies the Baco Island Shelf (BIS). Between the BIS and the Verde Island Shelf (VIS) lies the deeper Verde-Baco-Island Saddle (VBIS). The VBIS is ∼6-km long, ∼1.2-km wide, and has a maximum depth of less than 350 m. Both CS and VBIS connect the waters of Verde Island Passage and the Tayabas Bay, which lies east of the South Pass.

 4.3 The offshore trace of the Aglubang River Fault (ARF)

The previously described SW-facing, steep (10–45°) linear slope correlates with one of the nodal planes (Strike: N 21 W; Dip: 70) defined by the focal mechanism solution of the 15 November 1994 Mw 7.1 earthquake (Dziewonski et al., 1981; Ekstrom et al., 2012). This earthquake, associated with the Aglubang River Fault, resulted in a ∼35-km total onshore surface rupture in Oriental Mindoro (PHIVOLCS, 1994). As such, the NW-trending trough and scarp are interpreted as the approximate fault line scarp of the offshore extension of the Aglubang River Fault (ARF). This interpretation is in agreement with the mapped offshore extension of the ARF by Sarmiento and others (2022).

 4.4 The Malaylay Submarine Canyon System (MSCS)

A prominent feature identified within the South Pass is the NW-trending, 8-km long Malaylay Submarine Canyon System (MSCS; Figure 5), eponymously named due to its proximity to Malaylay Island in Oriental Mindoro. The MSCS can be classified as a headless, shelf-incising canyon due to the absence of a shelf valley connecting it to the adjacent outlets of the Baco and Malaylay rivers. The eastern side of the MSCS is structurally controlled by the Aglubang River Fault (ARF) fault line scarp described in the previous section.

Three-dimensional view of the offshore Aglubang River Fault and the Malaylay Submarine Canyon System (MSCS). The red dots show the location of depth comparisons with the thickness of eroded material in meters. KnC = Kanan Submarine Canyon; KlC = Kaliwa Submarine Canyon.

The MSCS is divided into a lower and upper course. The ∼4.5-km wide upper course is composed of two canyons, namely, the Kaliwa (KlC) and Kanan (KnC) canyons. These canyons are respectively equivalent to the Malaylay and Baco canyons described by Sequeiros and others (2019), with the names modified to avoid confusion with the MSCS.

Both the Kaliwa and Kanan canyons exit into the ∼6-km long, ∼1-km wide lower course dubbed as the Main Canyon by Sequeiros and others (2019). The Main Canyon eventually exits into the Verde Passage Basin at depths of ∼350 m.

 4.5 Submarine mass failures in the South Pass

A total of 15 probable submarine mass failures (SMFs) were identified within the South Pass, based on the occurrence of bight-shaped headwall scarps (Figure 4). Their corresponding head scarp lengths, areas, transport directions, and relative ages are also listed in Table 4. These SMFs were identified along the shelf edges of the eastern Northern Mindoro Shelf (eNMS), and within the Malaylay Submarine Canyon System (MSCS), the Verde Island Shelf (VIS), the Baco Island Shelf (BIS), the Calapan Saddle (CS), and the Verde–Baco Islands Saddle (VBIS).

The 15 mapped probable SMFs in the South Pass.

SMF ID Location a Latitude Longitude Length (m) Area Transport direction Relative age (post-1950s?) b
1 eNMS 13.43 121.17 880.72 0.103 N 09 E No
2 eNMS 13.42 121.16 910.32 0.115 N 24 W No
3 eNMS 13.42 121.15 636.40 0.053 Due N No
4 eNMS 13.42 121.14 413.25 0.034 N 8 W No
5 MSCS 13.42 121.12 2092.89 1.430 Due N Yes (9; 24–105 m)
6 BIS 13.47 121.16 603.36 0.056 N 07 E Yes (1; 31 m)
7 VIS 13.54 121.17 4,065.74 1.969 S 78 E N/A
8 VIS 13.55 121.13 2,702.10 0.822 N 44 E N/A
9 VIS 13.58 121.14 1,490.99 0.427 N 23 W N/A
10 VIS 13.56 121.03 1,244.12 0.408 S 55 W N/A
11 CS 13.43 121.13 432.05 0.041 N 63 E N/A
12 BIS/VBIS 13.51 121.18 409.66 0.036 N 13 W N/A
13 BIS/VBIS 13.51 121.18 415.55 0.036 Due N N/A
14 BIS/VBIS 13.51 121.18 500.78 0.049 N 56 W N/A
15 BIS 13.48 121.12 556.80 0.099 N 70 E N/A

eNMS–eastern Northern Mindoro Shelf; MSCS–Malaylay Submarine Canyon System; BIS–Baco Islands Shelf; VIS–Verde Island Shelf; CS–Calapan Saddle; VBIS–Verde–Baco Islands Saddle.

Based on comparison of 2011 bathymetry and 1950s point depth sounding data. If Yes, the number of comparison points and the range of net erosion values are also indicated. No means net erosion is <10 m. N/A means no comparison points were plotted downslope of the SMF.

The row in bold empahsizes the San Andres SMF, which is focused in this study.

The largest of the probable SMFs in terms of area are SMF 8 (0.822 km2), SMF 5 (1.430 km2), and SMF 7 (1.969 km2). Of the three, only SMF 5 can be reliably considered as having occurred after 1950s, with net erosion ranging from 24 to 105 m, based on nine measurements (Figures 5, 6). SMF 6 (0.056 km2) also shows significant net erosion of 31 m, although based only on one measurement. SMFs 1 to 4 show insignificant net erosion (<10 m), which may indicate their relatively older age. On the other hand, no comparison points were available for SMFs 7 to 15.

Initial sea surface displacement, tsunami propagation, and maximum wave height distribution of the EQ-only model. The Aglubang River Fault is shown in red, segmented into two sub-faults as described in Table 5.

Aside from its size and availability of net erosion measurements, SMF 5 is particularly notable due to its proximity to Barangay San Andres in Baco, Oriental Mindoro, which received the most damage from the tsunami. We therefore focus the following discussions on SMF 5, which we formally refer to as the San Andres SMF (SASMF).

 4.6 The San Andres submarine mass failure (SASMF) and the 1994 Mindoro tsunami

The proposed occurrence of the San Andres SMF (SASMF) has been based on the net erosion derived from the comparison between two, widely spaced bathymetry datasets. The question, therefore, is whether the observed net erosion was a product of a single, catastrophic event (e.g., submarine mass failure) or a series of smaller and/or gradual processes (e.g., turbidity currents).

The occurrence of erosional events within the MSCS since 1997 has been documented and studied by Sequeiros and others (2019), in relation to the destructive effects of erosion on a gas pipeline installed in 2001. Throughout a 21-year duration (January 1997 to September 2018) with a total of 27 pipe and aerial surveys, only two significant erosional events were documented within the MSCS, both of which were attributed to typhoon-triggered turbidity currents. Of the 11 typhoons that have passed over or within the vicinity of the MSCS, only typhoons Reming (Durian) in 2006 and Nina (Nock-Ten) in 2016 caused turbidity currents, resulting in pipe displacements and damage to protective rock berms. Seabed erosion due to these turbidity currents occur within established gullies and canyons of the MSCS. Erosion thickness generally does not exceed 5 m, except in some localized zones where it may reach up to 10 m. Moreover, the likelihood of turbidity currents within the MSCS due to hyperpycnal outflows from the Baco and Malaylay rivers has been ruled out, owing to insufficient sediment concentration even during annual monsoon floods.

We can then reasonably constrain the occurrence of the SASMF between the 1950s and 1997. Also, by extrapolating the frequency of typhoon-triggered turbidity currents within the MSCS as discussed above (∼1 in 11.5 yrs), we estimate four to five erosional events between the 1950s and 1997. Assuming a 5-m average thickness of eroded material per event, the cumulative erosion thickness is roughly 20–25 m. This estimate is on the lower end of the net erosion observed for the SASMF, which ranges from 24–105 m. Note that the simple estimate presented above assumes a negligible sediment input from the Baco and Malaylay rivers, which is unlikely as sediment yields of 0.336 and 0.351 MT/yr are estimated for both rivers, respectively (Sequeiros et al., 2019). Furthermore, in contrast to the much wider extent of the SASMF (∼1.43 km2), the extent of the turbidity currents’ influence is restricted within existing gullies and canyons within the MSCS.

While further constraints on the timing of the SASMF are unavailable at present, we hypothesize that the SASMF was likely triggered by the Mw 7.1 1994 Mindoro earthquake. This idea is premised on the following: 1) the 1994 Mindoro earthquake is within the possible window of occurrence for the SASMF, as already established above; and 2) modeling the 1994 Mindoro tsunami with the SASMF as an additional source mechanism approximates the observed wave height distribution better than the model based on the 1994 earthquake alone. We further elaborate on the second point in the following sections.

 5 Numerical modeling of the 1994 Mindoro tsunami 5.1 Modeling methods

The two models (EQ-only and EQ+SMF) for the 1994 Mindoro tsunami were calculated using JAGURS (Baba et al., 2015). JAGURS is a numerical code that calculates the propagation and inundation of tsunami waves over gridded terrain, based on the linear and nonlinear Shallow Water Equations (SWE), and Boussinesq Equations to account for wave dispersion (Baba et al., 2015).

5.1.1 Preparation of terrain data

The 20-m digital bathymetry and the 5-m IFSAR topography data from NAMRIA were used as terrain inputs for both models. Merging of the two datasets and interpolation for areas with unavailable data were done in QGIS3, using the tension spline interpolation method in the SAGA raster terrain module Close Gaps. Final preparation of the terrain data for JAGURS was done using GDAL and GMT.

 5.1.2 EQ-only model source mechanism

The source mechanism for the EQ-only model was based on the 15 November 1994 Mw 7.1 Mindoro earthquake along the right-lateral Aglubang River Fault (ARF). The initial sea surface displacement was calculated in JAGURS, based on the surface deformation due to a finite rectangular source (Okada, 1985). The effects of horizontal displacement (Tanioka and Satake, 1996) were also included in the calculations.

The fault parameters used in the calculations are shown in Table 5. We decided to segment the ARF into two sub-faults to reflect both the mapped offshore fault line scarp in this study and the mapped onshore fault rupture by PHIVOLCS (1994). As such, the strike of the sub-faults followed the general trends of the offshore fault line scarp and onshore fault rupture, while the downdip width is calculated using the scaling relationship equation based on moment magnitude by Wells and Coppersmith (1994). The dip angle (70°) is based on the GCMT focal mechanism of the Mw 7.1 earthquake (Dziewonski et al., 1981; Ekstrom et al., 2012). This dip angle generally coincides with the mapped offshore fault line scarp when projected from the hypocentral location (13.5 N, 121.1 E, 7–12 km depth) of the earthquake as reported by PHIVOLCS (1994). The slip amount was based on the maximum vertical displacement (1.9m) observed onshore at Baruyan Lake (PHIVOLCS, 1994).

Fault parameters for the Aglubang River Fault as used for the earthquake source mechanism. The two sub-faults are based on the geomorphological expression of the offshore extension (this study), and the mapped onshore rupture (PHIVOLCS, 1994) as labeled in Figures 6, 7.

Id Latitude a Longitude a Depth a , c (km) Strike Dip d Length (km) Width e (km) Rake d Slip b (m)
1 13.416 b , e 121.136 b , e 0 330 e 70 15.6 e 14 -178 2
2 13.245 b 121.179 b 0 342 b 70 19.2 b 14 -178 2

Location of the upper right corner of the rectangular fault plane.

Field and instrumental data from PHIVOLCS (1994).

Value selected as 0 for the surface-rupturing earthquake.

Global Centroid Moment Tensor (GCMT) focal mechanism solution (Dziewonski et al., 1981; Ekstrom et al., 2012).

Scaling relationships by Wells and Coppersmith (1994).

Derived from submarine geomorphological mapping (this study).

 5.1.3 EQ+SMF model source mechanism

A combined source mechanism is modeled based on the 1994 earthquake and SMF identified from submarine geomorphological mapping. For the earthquake component, the initial sea surface displacement from the EQ-only model was also used, as described in the previous section.

For the SMF component, we used the semi-empirical predictive equation by Sabeti and Heidarzadeh (2022), which was tested for the 1994 Skagway, Alaska and the 1998 Papua New Guinea submarine landslide-generated tsunamis. The equation relates the SMF volume (V) in km3 and the initial submergence depth (d) in m to the maximum initial wave amplitude (ηmax) in m, as follows: η max = 50.67 ( V d ) 0.34

The parameters used in the calculation of the SMF-generated wave are shown in Table 6. The initial sea surface displacement from the earthquake and SMF-generated wave were added to produce the EQ+SMF model initial conditions, based on the assumption that the tsunamigenic SMF was coseismic.

SMF parameters used in this study.

SMFa
Area (A) 1.43 sq. Km
Volume (V) 0.0483 cu. Km
Initial submergence depth (d) 10 m
SMF-generated tsunami wave
Amplitude (ƞ max )b 8.27 m
Radius (r) a 500 m
Central location (latitude, longitude) a 13.42176 N, 121.12342 E

Derived from submarine geomorphological mapping (this study).

Calculated using the equation by Sabeti and Heidarzadeh (2022).

 5.1.4 Parameters for tsunami propagation and inundation modeling

The numerical models used a terrain grid resolution of 0.00015° (∼15 m), with model bounds (NW and SE corners) of 13.65015 N, 120.9000 E and 13.4000 N, 121.19970 E. Total model duration was set at 15 min, with a calculation time step of 0.05 s. Snapshots of the model results were taken every 20 s. The effects of Coriolis motion and seawater compressibility were assumed negligible. The models used the linear shallow water equations to calculate the tsunami propagation.

 5.1.5 Calculation of modeled wave height statistics and error

We obtained modeled wave height statistics (i.e., mean, median, minimum, maximum, first quartile, third quartile, standard deviation) for 22 shoreline sections (SS1—SS22) around northern Oriental Mindoro, Baco Islands, and Verde Island.

We also compared the modeled wave heights (Wm) with 11 observed wave height measurements (Wo; see Supplementary Table S1) by calculating the root mean squared error (RMSE). The root mean squared error (RMSE) for both models was calculated using the equation, where the number of observations (n) is 11: R M S E = Σ ( W m − W o ) 2 n

5.2 The EQ-only and EQ+SMF tsunami models 5.2.1 Tsunami propagation

As previously described, we modeled two scenarios (EQ-only and EQ+SMF) for the 1994 Mindoro tsunami. We used the San Andres SMF as basis for the SMF component of the EQ+SMF model. The initial sea surface displacement, snapshots of tsunami propagation, and maximum wave height distributions for the EQ-only and EQ+SMF models are shown in Figures 6, 7, respectively.

Initial sea surface displacement, tsunami propagation, and maximum wave height distribution of the EQ+SMF model. The Aglubang River Fault is shown in red, segmented into two sub-faults as described in Table 5. The amplitudes of the SMF-generated wave crests and troughs are also labeled accordingly. In the last panel, note the focused distribution of extreme wave heights (>4 m) at Baco and Baco Islands, which is consistent with observed data, and is absent in the EQ-only model.

For the EQ-only model, initial sea surface displacements (Figure 6) ranged from -0.57 to 0.50 m. The sea surface displacements coincide with the locations of the insular slopes surrounding Verde Island, Baco Islands, and Calapan, which moved coseismically with the 1994 earthquake.

According to the EQ-only model, the movement of the ARF fault line scarp generated tsunami waves to the southwest, eventually reaching the coasts of Baco in 1 min; Puerto Galera in 2 min; and San Teodoro in 3 min. Meanwhile, Verde Island, Baco Islands, and Calapan were hit by tsunami waves within the first minute after the earthquake. Tsunami wave refraction around the northwestern point of Verde Island can also observed prominently 4 min after the earthquake. At least two distinct tsunami waves hit the coasts of Oriental Mindoro (Supplementary Figure S4), consistent with anecdotes from post-tsunami field surveys (PHIVOLCS, 1994; Imamura et al., 1995).

On the other hand, the EQ+SMF model is dominated by the SMF-generated wave, which is expected given its considerably larger amplitude (±8.27 m) compared to the earthquake-generated waves (Figure 7). The SMF produced a leading-depression wave towards Baco, consistent with eyewitness accounts (PHIVOLCS, 1994; Imamura et al., 1995). The SMF-generated wave spread out in all directions, with significant refraction towards the Baco Islands as it hit the Baco Island Shelf to the northeast. Focusing of tsunami waves were also observed leeward of the largest of the Baco Islands, as was similarly inferred from post-event runup measurments (Imamura et al., 1995).

Arrival times within the South Pass for both models are very short (<1–3 min), due to the near-field nature of the tsunami sources considered. However, the EQ+SMF model also shows farther propagation of tsunami waves north of Verde Island and towards Southern Luzon, as shown by the maximum wave height distributions (Figures 6, 7).

 5.2.2 Wave heights at shore

Figure 8 shows the comparison between modeled and observed wave heights at shore for 22 selected shoreline sections. The wave height statistics from the EQ-only and EQ+SMF models for the 22 selected shoreline sections are also listed in Supplementary Table S2.

Comparison between observed and modeled (EQ-only and EQ+SMF) wave heights. The upper panels show the locations of the shoreline sections, while the lower panels show the box plots of modeled wave heights for the EQ-only and EQ+SMF tsunami models. Observed wave heights are derived from run-up measurements by Imamura and others (1995), based on the wave-height-at-shore and run-up empirical relationship equation by Smart and others (2016).

For the EQ-only model, the predicted median wave heights for San Andres (SS9) and Baco Islands (SS20) are 0.63 m and 0.70 m, respectively. On the other hand, the predicted median wave heights from the EQ+SMF model are 9.78 m for San Andres, and 2.00 m for Baco Islands.

Furthermore, calculated RMSE values (Table 7) are consistently lower for the EQ+SMF model versus the EQ-only model. Using median wave heights at shore values, the RMSE for the EQ-only and EQ+SMF models are 3.21 m and 1.60 m, respectively. Using mean wave heights at shore values, the RMSE for the EQ-only and EQ+SMF models are 3.19 m and 1.72 m, respectively. Finally, using all wave heights at shore values, the RMSE for the EQ-only and EQ+SMF models are 3.20 m and 2.34 m, respectively.

Root mean squared error (RMSE) calculations for EQ-only and EQ+SMF models using mean, median, and all predicted wave heights at shore values.

( W m − W o ) 2 RMSE
n SS1 SS4 SS6 SS9 SS10 SS12 SS16 SS17 SS18 SS20 SS21
EQ-only (mean) 11 1.38 0.69 0.34 48.84 12.65 5.54 3.40 3.92 0.95 21.13 12.99 3.19
EQ-only (median) 11 1.54 0.76 0.36 48.81 12.63 6.24 3.47 3.91 0.75 21.41 13.24 3.21
EQ-only (all) 341 3.20
EQ+SMF (mean) 11 0.04 0.22 0.01 13.90 0.36 1.74 1.83 1.38 0.20 9.13 3.73 1.72
EQ+SMF (median) 11 0.00 0.21 0.04 4.70 0.75 2.51 2.27 1.51 0.26 11.00 4.85 1.60
EQ+SMF (all) 341 2.34
5.3 The M<sub>w</sub> 7.1 earthquake alone is not enough to produce the 1994 Mindoro tsunami

The results of our EQ-only model produced higher wave heights than those predicted by the earlier model of Tanioka and Satake (1996). Their model predicted wave heights of 0.07 m, 0.09 m, and 0.09 m for Baco, Wawa, and Calapan, respectively. Meanwhile, our EQ-only model predicted median wave heights at shore of 0.63 m (SS9), 0.74 m (SS12), and 0.98 m (SS16), respectively. These higher wave heights are likely due to the use of bathymetry data with a higher resolution. Lower resolutions smooth out bathymetric features such as insular shelves and slopes. This likely results in lower predicted wave amplitudes due to horizontal displacement of the seafloor, which is a significant contributor to initial tsunami wave heights. On the other hand, both models also predict higher wave heights for Calapan than in Baco, which is contrary to observed wave heights (2.83 m and 7.61 m, respectively).

More importantly, the EQ-only model severely underestimates the observed wave heights at shore for the 1994 Mindoro tsunami. For the EQ-only model, the predicted median wave heights for San Andres (SS9) and Baco Islands (SS20) are 0.63 m and 0.70 m, respectively. In contrast, observed wave heights at shore for San Andres and Baco Islands are 7.61 m and 5.32 m, respectively. Overall, the results of the EQ-only model clearly show that the Mw 7.1 1994 Mindoro earthquake alone is not enough to produce the wave height magnitudes and distribution observed for the 1994 Mindoro tsunami.

 5.4 The 1994 Mindoro tsunami has a combined earthquake and submarine mass failure source mechanism

The comparison between both modeled (EQ-only and EQ+SMF) and observed data wave height at shore distributions (Figure 8) shows that the EQ+SMF model approximates the observed data better than the EQ-only model. Quantitatively, this is also indicated by the lower RMSE values for the EQ+SMF model compared to the EQ-only model.

Considering the sample of the modeled wave height at shore measurements as a whole (n=341), the EQ-SMF model has a RMSE of 2.34 m versus the EQ-only model with a RMSE of 3.20 m. While this already indicates better predictions of the EQ+SMF model over the EQ-only model, the RMSE calculations for EQ+SMF are drastically affected by the large spread of predicted wave height at shore values specifically at SS9 (San Andres), with a standard deviation of 4.68 m. In contrast, the maximum spread of wave height at shore values for the EQ-only model is at SS12, with a standard deviation of 0.57 m. The larger spread observed for the EQ+SMF model is likely a function of the higher predicted wave heights compared to the EQ-only model. Moreover, constructive interference of waves would result in very high wave amplitudes at certain locations. The most extreme example is at SS9, which has predicted maximum and median wave heights at shore of 22.03 m and 9.78 m, respectively.

Furthermore, the RMSE of the EQ+SMF model improves if only the median or mean wave heights at shore values are considered in the calculation. RMSE values considering median and mean wave heights are 3.21 m and 3.19 m, and 1.60 m and 1.72 m, respectively. The two-fold increase in RMSE values between the EQ+SMF and EQ-only models clearly supports our hypothesis that the 1994 Mindoro tsunami is an earthquake-triggered, submarine mass failure tsunami.

 6 Conclusion

In this study, we combined submarine geomorphological mapping and numerical modeling to constrain and describe a more probable mechanism for the 1994 Mindoro tsunami. We have identified and mapped a previously undescribed submarine mass failure within the Malaylay Submarine Canyon System (MSCS), which we refer to as the San Andres submarine mass failure (SASMF). Furthermore, we have established evidence that constrains its occurrence between the 1950s and 1997, and ruled out turbidity currents as triggering mechanisms for extensive submarine mass failure. We therefore infer that the San Andres SMF is concurrent with the Mw 7.1 1994 Mindoro earthquake, and as such, is an important contributing mechanism to the resulting 1994 Mindoro tsunami.

To test this hypothesis, we modeled two scenarios for the 1994 Mindoro tsunami; the first was based on the Mw 7.1 quake alone (EQ-only model), while the second one included a submarine mass failure component based on the San Andres SMF (EQ+SMF model). The results of the numerical modeling show that the EQ+SMF model predicted wave height at shore values closer to observed data, compared to those predicted by the EQ-only model. As such, the results of our numerical modeling support the likelihood of a combined earthquake and submarine mass failure mechanism for the 1994 Mindoro tsunami.

On the other hand, significant improvements on the EQ+SMF model can and should still be made, as evident from the RMSE values, which are still in the same order of magnitude as the wave height measurements. Such discrepancies are expected given the limitations of how we calculated the submarine mass failure (SMF) component of the source mechanism. The calculations for the SMF component, based only on the volume and initial submergence depth, do not take into account the transport mechanics (i.e., translational, slumping, debris flow, turbidity flow) and dynamic aspects (i.e., velocity, acceleration) of the submarine mass failure (Harbitz et al., 2006; Masson, et al., 2006; Løvholt et al., 2015). Constraining these parameters requires a more detailed characterization of the San Andres SMF, which is beyond the scope of this work. Moreover, because bathymetry data obtained close to the 1994 tsunami is unavailable, we have opted to use the 2011 bathymetry data in the numerical modeling. This simplification has introduced another source of error, as changes in the bathymetry due to erosional and depositional processes would have definitely occurred between 1994 and 2011.

While past research on Philippine tsunamis have been focused on earthquake-generated mechanisms, our findings highlight the importance of including SMF-generated tsunamis in assessing the country’s coastal hazards. With the increasing availability of high-resolution bathymetry data not just for the inland and territorial seas of the Philippines but also globally, our methods on submarine geomorphological mapping and numerical modeling can be replicated to map and identify past and potential SMFs and estimate the magnitude and extent of the resulting tsunamis. This is particularly important for SMFs forming along shelves proximal to populated coastal areas like the San Andres SMF, which render conventional tsunami warning systems insufficient due to very short tsunami arrival times.

 Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

 Author contributions

AR conducted the submarine geomorphological mapping and numerical modeling, and wrote the paper. NR provided supervision in the conceptualization of the study, data analysis, and interpretation. KS, NR, TB, LN, RM, and IN contributed insights on the earthquake and submarine mass failure source parameters. KS, TB, and NC provided technical training and feedback on the numerical modeling of the 1994 Mindoro tsunami. NR and AR obtained funding for the study. NR, KS, TB, LN, IN, and RM were involved in the manuscript revisions.

 Funding

This study was primarily funded by the Department of Science and Technology - Philippine Council for Industry, Energy and Emerging Technology Research and Development (DOST-PCIEERD) Project No. 4233, co-implemented by the University of the Philippines–National Institute of Geologic Sciences (UP-NIGS) and Philippine Institute of Volcanology and Seismology (PHIVOLCS). Additional support was also provided by UP through NIGS Research Grants 2020 and 2022.

The authors would like to thank the National Mapping Resource and Information Authority (NAMRIA) for providing the bathymetry and topography datasets; to the local government units of Oriental Mindoro, including the municipalities of Baco, Calapan, San Teodoro, and Puerto Galera, for assistance during our fieldwork in February 2019; to external reviewers ABSR and LGR; and to the Geomorphology and Active Tectonics Research (GEAR) Laboratory. Gratitude is also extended to the contributors behind JAGURS, QGIS, SAGA, GDAL, and GMT.

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.1067002/full#supplementary-material

References Baba T. Takahashi N. Kaneda Y. Ando K. Matsuoka D. Kato T. (2015). Parallel implementation of dispersive tsunami wave modeling with a nesting algorithm for the 2011 Tohoku Tsunami. Pure Appl. Geophys. 172 (12), 3455–3472. 10.1007/s00024-015-1049-2 Badillo V. L. Astilla Z. C. (1978). Moro Gulf tsunami of 17 August 1976. Quezon City: Manila Observatory. Bautista M. L. P. Bautista B. C. Salcedo J. Narag I. (2012). Philippine tsunamis and seiches (1589 – 2012): A catalogue of tsunami events in the Philippines. Quezon City: Philippine Institute of Volcanology and Seismology (PHIVOLCS) 978-971-92395-1-2. Chiba T. Kaneta S. Suzuki Y. (2008). Red relief image map: New visualization method for three dimensional data. Remote Sens. Spatial Inf. Sci. 37. Claro S. M. D. Ramos N. T. Fernando A. G. S. Ishimura D. Switzer A. D. (2021). Sedimentological evidence of washover deposits from extreme wave events in Zamboanga del Sur, Mindanao, southern Philippines. Mar. Geol. 438, 106535. 10.1016/j.margeo.2021.106535 Conrad O. Bechtel B. Bock M. Dietrich H. Fischer E. Gerlitz L. (2015). System for automated geoscientific analyses (SAGA) v. 2.1.4. Geosci. Model Dev. 8, 1991–2007. 10.5194/gmd-8-1991-2015 De Reu J. Bourgeois J. Bats M. Zwertvaegher A. Gelorini V. De Smedt P. (2013). Application of the topographic position index to heterogeneous landscapes. Geomorphology 186, 39–49. 10.1016/j.geomorph.2012.12.015 Dziewonski A. M. Chou T. A. Woodhouse J. H. (1981). Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res. 86, 2825–2852. 10.1029/JB086iB04p02825 Ekström G. Nettles M. Dziewoński A. M. (2012). The global CMT project 2004–2010: Centroid-moment tensors for 13, 017 earthquakes. Phys. Earth Planet. Interiors 200-201, 1–9. 10.1016/j.pepi.2012.04.002 Fryer G. J. Watts P. Pratson L. F. (2004). Source of the great tsunami of 1 April 1946: A landslide in the upper aleutian forearc. Mar. Geol. 203 (3-4), 201–218. 10.1016/s0025-3227(03)00305-0 GDAL/OGR contributors (2020). GDAL/OGR geospatial data abstraction software library. Retrieved from: https://gdal.org . Gervasio F. C. (1971). Geotectonic development of the Philippines. J. Geol. Soc. Phil. 25 (1), 18. Guisan A. Weiss S. B. Weiss A. D. (1999). GLM versus CCA spatial modeling of plant species distribution. Plant Ecol. 143, 107–122. 10.1023/a:1009841519580 Gusman A. R. Supendi P. Nugraha A. D. Power W. Latief H. Sunendar H. (2019). Source model for the tsunami inside Palu Bay following the 2018 Palu earthquake, Indonesia. Geophys. Res. Lett. 46, 8721–8730. 10.1029/2019gl082717 Hamburger M. W. Cardwell R. K. Isacks B. L. (1983). Seismotectonics of the northern Philippine island arc. Am. Geophys. Union Monogr. 27, 1. Harbitz C. B. Løvholt F. Pedersen G. Masson D. G. (2006). Mechanisms of tsunami generation by submarine landslides: A short review. Nor. J. Geology/Norsk Geol. Forening 86 (3). Hayes D. E. Lewis S. D. (1984). A geophysical study of the Manila Trench, Luzon, Philippines 1. Crustal structure, gravity and regional tectonic evolution. J. Geophys. Res. 89, 9171–9195. 10.1029/jb089ib11p09171 Heidarzadeh M. Muhari A. Wijanarto A. B. (2018). Insights on the source of the 28 September 2018 Sulawesi Tsunami, Indonesia based on spectral analyses and numerical simulations. Pure Appl. Geophys. 176 (1), 25–43. 10.1007/s00024-018-2065-9 Heidarzadeh M. Satake K. (2017). Possible dual earthquake–landslide source of the 13 November 2016 Kaikoura, New Zealand Tsunami. Pure Appl. Geophys. 174 (10), 3737–3749. 10.1007/s00024-017-1637-4 Imamura F. Synolakis C. E. Gica E. Titov V. Listanco E. Lee H. J. (1995). Field survey of the 1994 Mindoro island, Philippines tsunami. Pure Appl. Geophys. 144 (3-4), 875–890. 10.1007/bf00874399 Karig D. E. (1983). Accreted terranes in the northern part of the Philippine archipelago. Tectonics 2, 211–236. 10.1029/tc002i002p00211 Lagmay A. M. Tejada M. L. Pena R. E. Aurelio M. A. Davy B. David S. (2009). New definition of philippine plate boundaries and implications to the philippine mobile belt. J. Geol. Soc. Philipp. 64, 17. Løvholt F. Pedersen G. Harbitz C. B. Glimsdal S. Kim J. (2015). On the characteristics of landslide tsunamis. Phil. Trans. R. Soc. A 373 (2053), 20140376. 10.1098/rsta.2014.0376 Masson D. G. Harbitz C. B. Wynn R. B. Pedersen G. Lovholt F. (2006). Submarine landslides: Processes, triggers and hazard prediction. Phil. Trans. R. Soc. A 364 (1845), 2009–2039. 10.1098/rsta.2006.1810 Mines and Geosciences Bureau (2010). Geology of the Philippines. Quezon city: Mines and Geosciences Bureau. Editors Aurelio M. A. Peña R. E. . 2nd Edition. Mitchell A. G. H. Hernandez F. de la Cruz A. P. (1986). Cenozoic evolution of the Philippine archipelago. J. Southeast Asian Earth Sci. 1, 3–22. 10.1016/0743-9547(86)90003-6 NCEI/WDS (2022). Global historical tsunami database. NOAA National Centers for Environmental Information. (Accessed on November 2022). 10.7289/V5PN93H7 Okada Y. (1985). Surface deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc. Am. 75 (4), 1135–1154. 10.1785/bssa0750041135 Okal E. A. Synolakis C. E. (2004). Source discriminants for near-field tsunamis. Geophys. J. Int. 158 (3), 899–912. 10.1111/j.1365-246x.2004.02347.x Ozawa A. Tagami T. Listanco E. L. Arpa C. B. Sudo M. (2004). Initiation and propagation of subduction along the Philippine Trench: Evidence from the temporal and spatial distribution of volcanoes. J. Asian Earth Sci. 23, 105–111. 10.1016/s1367-9120(03)00112-3 Pakoksung K. Suppasri A. Imamura F. Athanasius C. Omang A. Muhari A. (2019). Simulation of the submarine landslide tsunami on 28 september 2018 in Palu Bay, Sulawesi island, Indonesia, using a two-layer model. Pure Appl. Geophys. 176, 3323–3350. 10.1007/s00024-019-02235-y Pawluszek K. (2018). Landslide features identification and morphology investigation using high-resolution DEM derivatives. Nat. Hazards (Dordr). 96, 311–330. 10.1007/s11069-018-3543-1 PHIVOLCS (1994). 15 november 1994 Mindoro earthquake: Preliminary report of investigation. Quezon City, Philippines. QGIS (2022). QGIS geographic information system. QGIS Association. http://www.qgis.org. Rangin C. Spakman W. Pubellier M. Bijwaard H. (1999). Tomographic and geological constraints on subduction along the eastern Sundaland margin (South-East Asia). Bull. Geol. Soc. France 170, 775. Rangin C. Stephan J. F. Blanchet R. Baladad D. Bouysee P. Chen M. P. (1988). Seabeam survey at the southern end of the Manila trench. Transition between subduction and collision processes, offshore Mindoro Island, Philippines. Tectonophysics 146, 261–278. Issues 1–4. 10.1016/0040-1951(88)90095-9 Rangin C. (1991). The philippine mobile belt: A complex plate boundary. J. Southeast Asian Earth Sci. 6 (3-4), 209–220. 10.1016/0743-9547(91)90068-9 Sabeti R. Heidarzadeh M. (2022). A new empirical equation for predicting the maximum initial amplitude of submarine landslide-generated waves. Landslides 19, 491–503. 10.1007/s10346-021-01747-w Sarmiento K. J. S. Aurelio M. A. Flores P. C. M. Carrillo A. D. V. Marfito B. J. Abigania M. I. T. (2022). Seafloor structures and static stress changes associated with two recent earthquakes in offshore southern Batangas, Philippines. Front. Earth Sci. 9. 10.3389/feart.2021.801670 Sassa S. Takagawa T. (2018). Liquefied gravity flow-induced tsunami: First evidence and comparison from the 2018 Indonesia Sulawesi earthquake and tsunami disasters. Landslides 16, 195–200. 10.1007/s10346-018-1114-x Satake K. Tanioka Y. (2003). The july 1998 Papua New Guinea earthquake: Mechanism and quantification of unusual tsunami generation. Pure Appl. Geophys. 160 (10-11), 2087–2118. 10.1007/s00024-003-2421-1 Scarselli N. (2020). Submarine landslides – architecture, controlling factors and environments. A summary. Regional Geol. Tect. Princ. Geol. Analysis, 417–439. 10.1016/b978-0-444-64134-2.00015-8 Sequeiros O. E. Bolla Pittaluga M. Frascati A. Pirmez C. Masson D. G. Weaver P. (2019). How typhoons trigger turbidity currents in submarine canyons. Sci. Rep. 9 (1), 9220. 10.1038/s41598-019-45615-z Smart G. M. Crowley K. H. M. Lane E. M. (2015). Estimating tsunami run-up. Nat. Hazards (Dordr). 80 (3), 1933–1947. 10.1007/s11069-015-2052-8 Synolakis C. E. Bardet J. P. Borrero J. C. Davies H. L. Okal E. A. Silver E. A. (2002). The slump origin of the 1998 Papua New Guinea Tsunami. Proc. R. Soc. Lond. A 458 (2020), 763–789. 10.1098/rspa.2001.0915 Tanioka Y. Satake K. (1996). Tsunami generation by horizontal displacement of ocean bottom. Geophys. Res. Lett. 23 (8), 861–864. 10.1029/96gl00736 Van Den Eeckhaut M. Kerle N. Poesen J. Hervás J. (2012). Object-oriented identification of forested landslides with derivatives of single pulse LiDAR data. Geomorphology 173-174, 30–42. 10.1016/j.geomorph.2012.05.024 Von Huene R. Kirby S. Miller J. Dartnell P. (2014). The destructive 1946 Unimak near-field tsunami: New evidence for a submarine slide source from reprocessed marine geophysical data. Geophys. Res. Lett. 41 (19), 6811–6818. 10.1002/2014gl061759 Weiss A. D. (2001). Topographic position and landforms analysis. San Diego, CA: Poster Presentation, ESRI Users Conference. Wells D. Coppersmith K. (1994). New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. Bull. Seismol. Soc. Am., 84, 974–1002. Wilson J. P. Gallant J. C. (2000). “Digital terrain analysis,” in Terrain analysis - principles and applications. New York: John Wiley & Sons, 1. Yokoyama R. Pike R. (2002). Visualizing topography by openness: A new application of image processing to digital elevation models. Photogrammetric Eng. Remote Sens. 68 (3). Yumul G. P. Dimalanta C. B. Maglambayan V. B. Marquez E. J. (2008). Tectonic setting of a composite terrane: A review of the philippine island arc system. Geosci. J. 12 (1), 7–17. 10.1007/s12303-008-0002-0 Yumul G. P. Dimalanta C. B. Marquez E. J. Queaño K. L. (2009). Onland signatures of the palawan microcontinental block and philippine mobile belt collision and crustal growth process: A review. J. Asian Earth Sci. 34 (5), 610–623. 10.1016/j.jseaes.2008.10.002