Seismic and tsunami hazard potential of the Negros–Sulu megathrust, Philippines

The Negros–Sulu megathrust poses an imminent threat to coastal communities surrounding the Sulu Sea, but with limited information on past tsunamis, megathrust geometry, and locked segments, robust seismic and tsunami hazard assessments are hindered. To identify highly exposed coastal areas amid this knowledge gap, we estimated the Negros–Sulu megathrust source parameters using a structural-based segmentation and scaling-relation approach for tsunami modeling. A total of 18 sets of source parameters from six segments and three dip angle scenarios were considered with moment magnitude and average slip of M_w 8.0–8.9 and 1.68–6.23 m, respectively. Tsunami simulations were modeled in JAGURS, accounting for nonlinear shallow water equations, horizontal and vertical seafloor displacements, and Boussinesq dispersion effects. Coastal regions directly facing the segments have the highest exposures with <2 min arrival times, highlighting the major control of wave directivity and the need for rapid evacuation strategies. The Negros Trench generates up to 6 m wave height and 7 m s⁻¹ flow velocity, while the Sulu Trench produces up to 8 m and 6 m s⁻¹. Coastal areas with ≥2 m wave heights typically exhibit a concave morphology with a nearshore width interquartile range of 2–4.5 km. At wider nearshore widths (>20 km), wave dissipation results in lower wave heights (<2 m) that underscore the need for accurate nearshore bathymetry in tsunami modeling. As average slip is used in these models, we highlight that the maximum wave height and flow velocity are likely underestimated and recommend treating them as minimum values. Nonetheless, the relative exposure levels of coastal areas from these results provide a significant basis for regional-scale tsunami risk mitigation strategies and for searching tsunami deposits.

1 Introduction

Tsunamis brought about by megathrust ruptures are infrequent but have led to the most damaging coastal disasters in recorded history. The M_w 7.8–8.2 1923 Taisho Kanto (Satake, 2023), M_w 9.5 1960 Valvidia (Duke, 1960), M_w 9.2 1964 Great Alaska (Brocher et al., 2014), M_w 9.2 2004 Indian Ocean (Athukorala and Resosudarmo, 2005), and M_w 9.1 2011 Tohoku-Oki (Kajitani et al., 2013; Suppasri et al., 2012) events led to accumulated casualties of more than 380 000 and billions (USD) worth of damage. In the Philippines, Bautista et al. (2012) collated 41 positive tsunamigenic earthquake events from 1589 to 2012, primarily based on eyewitness accounts and previous reports. Of the 41 events (Fig. 1a), 20 can be directly linked to subduction zones, and three subduction-related events have a magnitude of at least 8.0: The M_s 8.0 1918 and M_w 8.1 1976 events along the Cotabato Trench with up to 8 and 9 m maximum tsunami wave height, respectively, and the M_s 8.2 1946 event along the Negros Trench but with a lower maximum tsunami wave of 2 m. The M_w 8.1 1976 is the most damaging historical tsunami event in the Philippines, characterized by a shallow rupture with a 25 km depth leading to estimated casualties of 3000–8000 (Acharya, 1978; Badillo and Astilla, 1978; Stratta et al., 1977). A washover deposit (766–675 cal BP) was also discovered along the coasts of Zamboanga del Sur facing the Cotabato Trench, suggesting an extreme wave event (tsunami or storm surge) prior to the 1976 earthquake and tsunami (Claro et al., 2021).

Figure 1 Philippine tectonic setting. (a) Estimated epicenters of positive tsunamigenic earthquakes and the associated eyewitness accounts of maximum wave heights in the affected coastal areas since 1589 (Bautista et al., 2012). The trace of the Philippine Fault Zone (PFZ) was adapted from PHIVOLCS. Depth contours of the Slab 2.0 model from Hayes et al. (2018). Arrows indicate the velocity vectors relative to no-net-rotation reference frame from inland GPS networks (Kreemer et al., 2014). Bathymetry and topography data from the 2024 General Bathymetric Chart of the Oceans (© GEBCO Compilation Group, 2024). (b) Crustal thickness (Laske et al., 2013), (c) sediment thickness (Straume et al., 2019), and (d) Bouguer anomaly (Balmino et al., 2012) across the archipelago.

One of the primary goals when inferring future megathrust ruptures is to establish seismic gaps, and coupling models with slip deficits are used to derive these locked segments from geodetic measurements (e.g., Chlieh et al., 2008; Métois et al., 2013, 2016; Michel et al., 2019; Moreno et al., 2010; Ozawa et al., 2012). However, with the current temporal limitation of geodetic records compared with the hundred-to-thousand-year recurrence of megathrust earthquakes, uncertainty remains regarding the likelihood and extent of ruptures. Even in subduction zones with low coupling, evidence of past ruptures and tsunamis was identified (e.g., Bécel et al., 2017; Debaecker et al., 2026; Witter et al., 2016), highlighting the limitations of solely relying on coupling models to infer future ruptures.

Coupling models have been constrained in the Manila, East Luzon, and Philippine Trench as well as along the Philippine Fault Zone (PFZ) (Galgana et al., 2007; Hsu et al., 2012, 2016) as GPS networks are dense in these regions. On the other hand, a sparse GPS network along the Negros–Sulu Trench inhibits identifying locked segments with slip deficits that may rupture. The current Slab 2.0 of global subduction zone geometry (Hayes et al., 2018) includes all trenches in the Philippines except the Negros–Sulu Trench, and while previous studies have provided tomography profiles of the Negros Trench (Fan et al., 2017; Fan and Zhao, 2019), the geometry of the Negros–Sulu megathrust remains uncertain.

Although earlier source estimation and tsunami modeling exist in this active margin (Nurashid et al., 2013; Okal et al., 2011; Pedersen et al., 2011; Salcedo and Hara, 2011; Zhao and Niu, 2022), the uncertainties of the slab geometry were not accounted for, as well as updated scaling equations more suitable for megathrust ruptures and hydrodynamic complexities such as nonlinear effects and Boussinesq dispersion. To fill these gaps and provide crucial information on the tsunami wave characteristics, this study estimated the megathrust earthquake source parameters in the Negros and Sulu subduction zones using a structural-based segmentation and scaling-relation approach. We produced three slab dip scenarios and used the estimated source parameters to model tsunamis and identify potentially exposed coastal areas across the coastal areas surrounding the Sulu Sea.

2 Seismotectonic Setting

The amalgamation of different terranes in the Philippine Mobile Belt is a consequence of the complex interaction among the Pacific, Philippine Sea, Eurasia, and Indo–Pacific plates, leading to its present configuration marked by oppositely dipping subduction zones and the PFZ across the archipelago (Aurelio, 2008; Aurelio et al., 2013; Hilde et al., 1977; Holloway and Hall, 1996; Yumul et al., 2003, 2008; Zahirovic et al., 2014). These subduction zones include the East Luzon Trough, Philippine, Manila, Negros–Sulu, and Cotabato trenches that can potentially generate great to epic (M>8) megathrust earthquakes (e.g., Li et al., 2016; Qiu et al., 2019; Ramos et al., 2017; Wu and Huang, 2009). Peñarubia et al. (2020) evaluated an empirical Gutenberg–Richter (GR) magnitude–frequency distribution (MFD) in these subduction zones. The b-values range from 0.691 to 1.724, based on the declustered seismicity catalog from the International Seismological Centre (ISC) and the Philippine Institute of Volcanology and Seismology (PHIVOLCS), where <1.0 indicates that large earthquakes contribute more to the long-term seismic moment rate (Nishikawa and Ide, 2014). In the study, the number of earthquakes in the declustered catalog in all trench segments was less than 300, with the majority of the segments having <50 (Peñarubia et al., 2020). In the Negros and Sulu subduction zones, only 7 and 21 earthquakes, respectively, were used to evaluate the GR MFD, with resulting b-values of 1.084 for Negros and 0.797 for Sulu. With longer and denser seismicity records, constraints on the GR MFD in this region will further improve. However, even with relatively denser declustered seismicity, rheological conditions and seismic behavior are not fully captured due to the complex factors affecting frictional heterogeneity (Chlieh et al., 2008; Fagereng and Beall, 2021; Fan and Zhao, 2021; Volpe et al., 2023; Yabe and Ide, 2018) and the short span of observations relative to the apparent hundred-to-thousand-year recurrence of megathrust earthquakes.

Inland GPS networks in the Philippines (Galgana et al., 2007; Hsu et al., 2016; Kreemer et al., 2014; Rangin et al., 1999) are mostly distributed in Luzon and eastern Philippines, revealing a west-northwest motion (Fig. 1a), while sparse GPS data in the southwestern Philippines show a south–southwest direction. The oceanic plates subducting beneath the five subduction zones have a crustal thickness ranging from 7 to 15 km (Fig. 1b), where relatively thicker crusts (>15 km) are linked to the island arcs across the Philippine Archipelago (Laske et al., 2013). PHILCRUST 3.0 provides a higher resolution inland crustal thickness varied by different island growth rates and terranes (Parcutela et al., 2020). Based on the global dataset from Straume et al. (2019), the majority of marginal seas and active margins in the Philippines are characterized by ≥1 km thick overlying sediments (Fig. 1c), except in the Philippine Sea. The Sulu Sea is particularly marked by high sediment thickness of up to 6 km. Estimates on the gravity variations in the Philippines from Balmino et al. (2012) show distinct high Bouguer anomalies (>300 mGal) of the subducting oceanic basins where the trench marks the boundary of the active margin and island arcs with relatively lower (<300 mGal) Bouguer anomaly signatures (Fig. 1d). The NW Sulu sub-basin, with the thickest sediments, has a lower Bouguer anomaly than the SE Sulu sub-basin. Past GPS campaigns estimated convergence rates of <60 mm yr⁻¹ along the Negros–Sulu subduction zone (Rangin et al., 1999; Simons et al., 1999).

3 Methodology

3.1 Constraining Megathrust Fault Parameters

Simplified rectangular fault models were adapted for the megathrust with parameters length (km), width (km), strike (°), dip angle (°), depth from the seafloor (km), rake (°), and slip amount (m). Insights from previous seismic tomography profiles (Fan et al., 2017; Fan and Zhao, 2019; Sevilla, 2011) on the geometry of the subducting slab reveal tapering of high P-wave velocity in the northern Negros Trench segment west of Panay Island, which provides limiting constraints on the length of the northern Negros Trench segment. Steepening of the subducting slab was also observed in the seismic tomography profiles around 100 km depth, although this is well beyond global seismogenic zones. Best-fit estimates on the slab dip of the shallow portion of the subducting slab derived from the tomography profiles within 50 km depth range from 9 to 30°. Focal mechanism solutions (ISC data from 1977 to 2023) of reverse-oblique rake values (10 to 170°) (Bondár et al., 2015; Di Giacomo et al., 2015, 2018; Storchak et al., 2013, 2015) also show a wide distribution of dip angles ranging from 9 to 90° (Fig. 2). To account for this uncertainty of the slab dip, we used three dip angles for all source parameters at 10, 20, and 30°.

Figure 2 Reverse-oblique focal mechanism solutions (1977–2023) of earthquakes along the Negros–Sulu subduction zone collated by the International Seismological Centre (ISC) (Bondár et al., 2015; Di Giacomo et al., 2015, 2018; Storchak et al., 2013, 2015). Earthquakes in the cross-sectional profiles (NT–NT^′ and ST–ST^′) are within 80 km depth and 200 km distance from each side of the profile. The possible interface of the megathrust fault plane is marked in the light blue zone, and dashed lines indicate the upper and lower boundaries of the megathrust rupture.

We utilized a structural-based segmentation in which the megathrust fault length is derived from the orientation variations of the frontal wedge (i.e., Nawanao and Ramos, 2023), and the width and bottom depth are limited by the seismogenic zone. These led to four main segments of the Negros–Sulu subduction zone, namely, NT1, NT2, ST1, and ST2. The NT1 and NT2 are marked by north-northeast and north-northwest orientations, respectively, while the ST1 and ST2 by north–northeast and northeast, respectively (Fig. 3a). Due to the uncertainty of the northern extent of NT1 where P-wave velocity profiles of Fan et al. (2017) and Fan and Zhao (2019) indicate tapering of the slab northward, we included longer (NT1-1) and shorter (NT1-2) rupture lengths for NT1. For the Sulu Trench, an additional longer segment (ST) that combines ST1 and ST2 is also modeled due to their similar northeast trend. This segmentation results in six subfaults, i.e., NT1-1, NT1-2, NT2, ST, ST1, ST2.

Figure 3 Megathrust fault parameterization of Negros–Sulu Trench segments. (a) A total of 18 sets of megathrust fault parameters were utilized for the tsunami modeling based on six segments and three dip angle scenarios. (b) Cross-section schematic diagram of the megathrust plane, where the width saturates at 196 km and is calculated using sine function. (c) Six cross-sectional segments of the Negros–Sulu subduction zone projected in the N–S profile. The colors of the megathrust plane follow those in Fig. 3a.

A shallow blind rupture is assumed for all scenarios with a top depth of 5 km and a rake of 90° (purely reverse slip). Based on the interquartile range (1.5 IQR) of 207 reverse-oblique earthquakes (Fig. 2), a threshold depth of 50 km was set as the bottom depth of the seismogenic zone. The width of the megathrust is calculated using sine function from the top and bottom depth, with the given the dip angles, and we adopted a saturation limit of 196 km as observed in global subduction zones (Allen and Hayes, 2017; Skarlatoudis et al., 2016; Tajima et al., 2013). The dip angles 10, 20, and 30° correspond to 196 (maximum width limit), 132, and 90 km, respectively (Fig. 3).

To account for the variability of length and width, the area (S) of the fault plane is used as input in the linear rupture-scaling equations (Eqs. 1–2) of Allen and Hayes (2017) to estimate the maximum M_w and the average (D_ave) slip (m). The linear area–M_w was preferred over the bilinear area–M_w scaling equation as the former produced slightly higher M_w, adopting a more conservative approach. A total of 18 sets of source parameters were constrained for elastic dislocation and tsunami modeling, with the ST (379 km) and NT1-1 (201 km) having the longest lengths (Table 1). Estimates of the moment magnitude from scaling equations range from M_w 8.0 to 8.9 (±0.27), where the ST segment with the longest fault length has the largest magnitude. These earthquake magnitudes correspond to an average slip of 1.68–6.23 m (±0.21).

$$\begin{array}{}\text{(1)}& {\displaystyle }{\log }_{\mathrm{10}}S=-\mathrm{3.63}+\mathrm{0.96}\times M_{\mathrm{w}}\text{(2)}& {\displaystyle }{\log }_{\mathrm{10}}{D}_{\mathrm{ave}}=-\mathrm{5.05}+\mathrm{0.66}\times M_{\mathrm{w}}\end{array}$$

Table 1 Negros–Sulu megathrust fault parameters of the six segments and three dip angle scenarios based on the scaling equations of Allen and Hayes (2017).

3.2 Tsunami Modeling

We utilized the 15 arcsec-resolution 2024 General Bathymetric Chart of the Oceans (GEBCO Compilation Group, 2024) bathymetry data as input for the tsunami simulation. To have a higher elevation accuracy of the topography, we resampled the interferometric synthetic aperture radar (IFSAR)-derived 5 m-resolution digital terrain model (DTM) (±1 m vertical accuracy) from the National Mapping and Resource Information Authority (NAMRIA) to 15 arcsec and merged it with the GEBCO bathymetry data. The resulting 15 arcsec IFSAR DTM–bathymetry raster served as the input depth in computing the initial seafloor and seawater displacement for tsunami modeling. As the entire coastal areas surrounding the Sulu Sea were considered for tsunami simulations, high-resolution subgrids for inundation models were not performed due to the absence of high-resolution nearshore bathymetry as well as restrictions in computational capacity and processing time.

The elastic deformation from the finite rectangular fault source is analytically solved through Green's function, assuming a homogeneous isotropic elastic half-space (Okada, 1985) and utilizing the average slip from the scaling equation. The tsunami wave propagation is modeled using two-dimensional nonlinear shallow water equations through the finite difference method (FDM) in a staggered grid deployed in the JAGURS numerical code (Baba et al., 2015, 2017). In JAGURS, both the initial vertical and horizontal seafloor displacement and Boussinesq dispersion effects are accounted for in the tsunami propagation. The modeling assumes a uniform surface roughness (Manning's coefficient = 0.025 s m${}^{-\mathrm{1}/\mathrm{3}}$), and we set a time step (Δt) of 0.5 s to maintain stability conditions based on the maximum depth (d), bathymetry resolution (Δx), and acceleration due to gravity (g) in Eq. (3) (Satake, 2002):

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }t\le {\displaystyle \frac{\mathrm{\Delta }x}{\sqrt{\mathrm{2}gd}}}\end{array}$$

The tsunami wave propagation time was set to 4 h, and 34 synthetic tide gauge stations (TGS) (Fig. 4 and Table S1 in the Supplement) were positioned across the coastal areas surrounding the Sulu Sea. The coastal areas surrounding the Negros–Sulu Trench are subdivided into 16 regions: I, Visayan Sea; II, Mindoro; III, Coron–Cuyo; IV, West Panay; V, Panay Gulf; VI, Sipalay–Bayawan; VII, Tañon–Cebu Strait; VIII, Bohol Sea; IX, Leyte Gulf; X, Palawan; XI, Cagayancillo; XII, Sindangan–Labason; XIII, Siocon–Sibuco; XIV, Basilan–Sulu; XV, Moro Gulf; and XVI, northeast Borneo. To represent the minimum arrival time, maximum wave height, and flow velocity across the coastlines, 6879 sampling points were set at 10 km intervals, each with a radius of 5 km.

Figure 4 The 16 coastal regions surrounding the Negros–Sulu subduction zone and the 34 synthetic tide gauge stations (TGS). (I) Visayan Sea; (II) Mindoro; (III) Coron–Cuyo; (IV) West Panay; (V) Panay Gulf; (VI) Sipalay–Bayawan; (VII) Tañon–Cebu Strait; (VIII) Bohol Sea; (IX) Leyte Gulf; (X) Palawan; (XI) Cagayancillo; (XII) Sindangan–Labason; (XIII) Siocon–Sibuco; (XIV) Basilan–Sulu; (XV) Moro Gulf; (XVI) northeast Borneo. Basemap from © GEBCO Compilation Group (2024).

4 Results

Tsunami Models

Figures 5 and S1–S5 (in the Supplement) show the representative wave propagation snapshots from rupture scenarios of six segments (NT1-1, NT1-2, NT2, ST, ST1, ST2) based on the 20° slab dip. The initial wave propagation from all ruptures is characterized by strong east–west and northwest–southeast wave crests in the Negros and Sulu Trench, respectively. Short-period wave dispersion and interferences are also prevalent as waves propagate toward shallower depths. In I, III, VII, VIII, and IX, the combined effects of dispersion, diffraction, and refraction caused significant deformation of wave crests as they propagated toward the coasts.

Figure 5 Tsunami wave (2 h) propagation from ST segment at 20° dip angle. Basemap from World Imagery Esri | Powered by Esri.

Coastal islands in direct proximity have <2 min arrival time, including V and VI for Negros Trench segments and XII, XIII, and XIV for Sulu Trench segments. Arrival times from Negros Trench segments in XI, XII, XIII, and northeast of X were about 10, 15, 30, and 45 min, respectively. As the initial wave approaches shallow depths, wave speed significantly drops as indicated by the compression of arrival time contours (Fig. 6). For instance, the first wave arrives in the southwest of X about 1 h after the initial displacement, while about 2 h in XVI. The effects of wave refraction, diffraction, and interference are also apparent in the deformed arrival time contours in the majority of coastal regions.

Figure 6 Tsunami wave arrival times from (a) NT1-1 and (b) ST segments. Contours have 5 min intervals. (I)–(XI) indicate the 16 coastal regions described in Fig. 4. Basemap from World Imagery Esri | Powered by Esri.

For ST and ST1 segments, the initial wave arrives in VI, XI, IV, and III by 10, 20, 25, and 45 min, respectively. For ST2 segments, arrival times are around 30 min in VI and XI, while around 35 min, 45 min, and 1 h in XVI, V, and III, respectively. Tsunami waves from Sulu Trench segments can also pass through the XIV and propagate across the Celebes Sea, arriving in XV at about 10 min (Fig. 6).

Waveforms from 34 synthetic tide gauge stations (Figs. 7–8 and S6–S7) also indicate the arrival as the first crest or trough. Relatively long-period waves are apparent in TGS 1 Coron, 2 Pandarochan Bay, 13 Himamaylan, 14 Bacolod, and 18 Cebu, while the majority are dominated by short-period waves due to nonlinear and dispersion effects. Tsunami waves from all segments did not propagate in TGS 22 Surigao City within the 4 h simulation. In TGS 33 Jolo, dissipation also led to no tsunami waveforms from NT and ST1, while the ST and ST2 ruptures caused coastal uplift.

Figure 7 Tsunami waveforms in TGS 1–16 from Sulu Trench segments (see Fig. 4 for locations).

Figure 8 Tsunami waveforms in TGS 17–34 from Sulu Trench segments (see Fig. 4 for locations).

The maximum wave height and flow velocity from six segments are shown in Fig. 9 (Figs. S8–S9) and Fig. 10 (Figs. S10–S11), respectively. The highest wave heights (>2 m) are observed in coastal areas parallel to the rupture segments, as also indicated in six tide gauge stations (13, 16, 28, 29, 30, 32). For NT1, tsunami waves reached up to 4–5 m in VI. For the NT2 segment, tsunami waves of up to 4 m are localized in VI. For ST and ST1 segments, the highest wave heights are located in XII (8 m) and XIII (5–7 m). Tsunami waves from ST segments propagate across the Celebes Sea, reaching 3–7 m in XV.

Figure 9 Maximum wave height (m) from the 4 h tsunami simulation of six segments with a megathrust dipping at 20°. Basemap from World Imagery Esri | Powered by Esri.

Figure 10 Maximum flow velocity (m s⁻¹) from the 4 h tsunami simulation of six segments with a megathrust dipping at 20°.

Coastal areas on the rear side and parallel to the trench segment also have relatively high tsunami waves. Negros Trench segments produced up to 2–3 m of tsunami wave heights in III and 1–2 m in XI and northeast of X. Sulu Trench segments also produced up to 3–4 m in X and up to 1–2 m in XI and III. The majority of the central Philippines (I, VII, VIII, IX), as well as in XVI, have wave heights of <1 m, except in some embayments and narrow channels where amplifications can reach 1–2 m.

The maximum tsunami wave height and flow velocity are coupled, i.e., coastal areas with the highest tsunami waves also have the highest flow velocity. For Negros Trench segments, the coastal areas with high tsunami waves in VI correspond to flow velocities of 3–7 m s⁻¹, while reaching up to 3 m s⁻¹ in III. For Sulu Trench segments, coastal areas in XII, XIII, XIV, and XV with the highest tsunami waves are exposed to a flow velocity of 3–6 m s⁻¹. In coastal areas of northeast Palawan, the flow velocity reached 3–4 m s⁻¹. The majority of coastal areas with <1 m tsunami wave height correspond to a maximum velocity of <1 m s⁻¹.

Table 2 summarizes the maximum, minimum, and average differences of the maximum tsunami wave height and flow velocity from megathrust segments at 10° dip compared to 20 and 30°. A shallow-dipping megathrust tends to produce the largest earthquakes (Bletery et al., 2016; Schellart and Rawlinson, 2013), hence the use of the 10°-dip rupture parameter as a reference. This tendency is highlighted by the positive average and maximum differences among parameters. However, in certain regions, steeper initial seawater displacement from rupture parameters with 20 and 30° dip angles produced waveforms that preferentially led to higher maximum wave height and flow velocity as represented by the negative minimum differences. The spatial distribution of the maximum wave height and wave velocity differences among parameters is indicated in Figs. 11 and S12–S14.

Figure 11 Maximum wave height differences (10–20 and 10–30°) comparing 10°-dip megathrust with those of 20 and 30° in the Sulu Trench (ST) segments.

Table 2 Statistics of the maximum wave height and flow velocity differences comparing 10°-dip megathrust with 20 and 30°.

For Negros Trench segments, the range of maximum difference is 0.934–2.507 m and 0.968–1.811 m s⁻¹ for maximum wave height and flow velocity, respectively, while the average differences are 0.008–0.104 m and 0.002–0.071 m s⁻¹. The minimum differences range from −0.529 to −1.141 m (maximum wave height) and −1.099 to −2.140 m s⁻¹ (maximum flow velocity). For Sulu Trench segments, the differences are higher than those from Negros Trench segments, with a range of 2.289–5.622 m (maximum wave height) and 2.369–4.509 m s⁻¹ (maximum flow velocity) for the maximum, 0.014–0.116 m and 0.003–0.017 m s⁻¹ for the average, and −0.481 to −1.890 m and −1.020 to −2.164 m s⁻¹ for the minimum.

5 Discussion

5.1 Influence of Wave Directivity, Coastal Morphology, and Nearshore Width on Maximum Wave Height

Wave directivity mainly controls the distribution of the highest tsunami wave heights, which are coastal regions directly facing (i.e., V, VI, XII, XIII, XIV, XV) and on the rear side (III, X) of the segments. As tsunami waves propagate in the interior central Philippines, wave interactions, refraction, diffraction, reflection, and nonlinear effects result in significant dissipation, with the majority of coastal areas in I, VII, XIII, and IX having <1 m wave height. However, the morphology of coastlines, especially along embayments, can amplify tsunamis due to constructive resonance, which can even generate up to 2 m in interior islands (Fig. 12b–f). Previous studies have also acknowledged the effects of morphology on tsunami wave resonance (e.g., Aranguiz et al., 2019; Roeber et al., 2010; Vela et al., 2014; Yamazaki and Cheung, 2011).

Figure 12 Representative sites highlighting the influence of nearshore bathymetry and coastal morphology on the maximum wave height among 18 segment scenarios. (a) Inset maps of (b)–(h). The white contours of the bathymetry from GEBCO have 50 m intervals. (b–d) southwest Negros (VI), (e) Labason (XII), (f) northern Bohol (VII), (g) Coron Islands (III), and (h) northeast Borneo (XVI). Basemap from World Imagery Esri | Powered by Esri.

The foremost islands can also function as a natural barrier for the rear islands as they diffract and attenuate the tsunami waves. This case is apparent for the higher tsunami waves in III compared to its rear islands (Fig. 12g). Another factor playing a role in the resulting tsunami waves is the nearshore width along the coasts, as nonlinear shallow-water waves dominate. Arbitrarily setting the distance between the −50 and 0 m depth contour as a reference for nearshore width in 6879 sampling points, the nearshore width interquartile range of the coastal areas with ≥2 m maximum wave height is 2–4.5 km with a maximum bound of 8 km (Fig. 13). On the other hand, coastal areas with <2 m maximum wave height have a nearshore width interquartile range of 2.8–10.6 km with a maximum bound of 22 km. All coastal areas with nearshore width ranging from 30 to 77 km also have <2 m maximum wave height. This indicates that as the nearshore width widens (e.g., >20 km), wave dissipation increases, reducing the tsunami wave as it reaches the shoreline. This is exemplified in the difference in distribution of maximum tsunami wave heights between the northern and southern coastline of XVI (Fig. 12h).

Figure 13 Maximum tsunami wave height from 6879 sampling points across the coastal regions in central and southern Philippines and the corresponding nearshore width (km) relative to 50 m depth. Boxplots show the interquartile range (1.5 IQR) of nearshore width between coastlines with ≥2 and <2 m maximum tsunami wave height. Orange circles indicate the outliers of the interquartile range.

5.2 Implications for Tsunami Hazard Potential

Accounting for all 18 megathrust rupture scenarios, Figs. 14, S15, and S16 present the maximum wave height (m), flow velocity (m s⁻¹), and minimum arrival times (min), respectively, in coastal areas surrounding the Sulu Sea. Of the 6879 sampling points, 8 % have a maximum wave height of ≥2 m. Sindangan–Labason (XII) and Siocon–Sibuco (XIII), Sipalay–Bayawan (VI), Moro Gulf (XV), and Basilan–Sulu (XIV) have the highest exposure of 5–8 m wave height and 5–7 m s⁻¹ flow velocity. Certain coastlines in Palawan (X), Panay Gulf (V), Coron–Cuyo (III), West Panay (IV), northeast Borneo (XVI), and Bohol Sea (III) are potentially exposed to 2–4 m wave height and 2–4 m s⁻¹ flow velocity. The majority of Tañon–Cebu Strait (VII), Mindoro (II), Leyte Gulf (IX), and Visayan Sea (I) coastal regions have a maximum wave height and flow velocity of <1 m and <1 m s⁻¹, respectively, with minor amplification of up to 2 m and 2 m s⁻¹ in a few areas. Comparisons of wave heights across 16 coastal regions are also shown in Fig. 15 with an average ranging from 0.1 to 2.1 m and standard deviations of ±0.06–0.77 m.

Figure 14 Map of coastal regions in central and southern Philippines showing the highest maximum wave height from the 18 megathrust rupture scenarios. The 6879 sampling points have a 10 km interval. Basemap from World Imagery Esri | Powered by Esri.

Figure 15 Comparisons of maximum wave heights among the 16 coastal regions from the models and eyewitness accounts (triangles and squares) of historical tsunamigenic earthquakes in the Sulu and Celebes Sea from Bautista et al. (2012). White circles and blue bars indicate the mean and maximum values, respectively, while the envelope indicates the 1σ and 2σ standard deviations. (I) Visayan Sea; (II) Mindoro; (III) Coron–Cuyo; (IV) West Panay; (V) Panay Gulf; (VI) Sipalay–Bayawan; (VII) Tanon–Cebu Strait; (VIII) Bohol Sea; (IX) Leyte Gulf; (X) Palawan; (XI) Cagayancillo; (XII) Sindangan–Labason; (XIII) Siocon–Sibuco; (XIV) Basilan–Sulu; (XV) Moro Gulf; (XVI) northeast Borneo. Rupture mechanisms: UP, upper plate; MT, megathrust; IP, intraplate; OR, outer rise. Subduction zones: NT, Negros Trench; ST, Sulu Trench; CT, Cotabato Trench.

The majority of coastal regions with the highest tsunami waves (XII, XIII, VI, XV, XIV) also have the shortest arrival times of <2 min due to their proximity to the rupture segments. This short arrival requires the need for quick evacuation strategies. In certain areas where established evacuation sites are far, vertical evacuation, i.e., reaching the closest elevated and structurally strong location, may be the most viable option. While the variability of coastal surface roughness is not accounted for in these models, pre-existing natural and man-made barriers may delay or inhibit inundation, providing additional time for evacuation.

We compared the maximum wave heights from the models to the eyewitness accounts of historical tsunamis in the Sulu and Celebes seas from Bautista et al. (2012). Although the source mechanisms of the majority of earthquakes were unknown, we can attribute them to megathrust, upper-plate, outer-rise, or intraplate ruptures, depending on the hypocenter when available and the position of the estimated epicenter relative to the trench, as labeled in Fig. 15. The estimated epicenter of the 1897 M_s 7.5 event was along the Sulu Trench with a 2 m maximum wave height in the majority of the Sulu Sea basin (III, IV, VI, X, XIII, XV, XVI) and up to 6 m in coastal region XIV. In III, IV, VI, X, XIII, XV, and XVI, the highest modeled tsunami wave in this study is slightly greater than the eyewitness accounts. However, in XIV where a 6 m wave height was estimated, the highest modeled tsunami wave is slightly lower at 5.4 m. Eyewitness accounts of the more localized tsunamis by the 2012 M_s 6.9, 1925 M_s 6.8, and 1948 M_s 8.2, reported up to 2 m wave height, which is lower than the highest modeled tsunami. In Moro Gulf (XV), six historical tsunamigenic earthquakes (2002 M_s 7.2, 1976 M_s 8.1, 1928 M_s 7.5, 1923 M_s 7.2, 1918 M_s 8, 1917 M_s 6.4) were generated along the Cotabato Trench (CT) with eyewitness accounts reporting tsunami wave heights ranging from 1 to 9 m. In this region, the highest modeled tsunami from the Sulu Trench reached 6 m. This indicates that this coastal region is not only susceptible to the Cotabato subduction zone rupture but also to that of the Sulu subduction zone. In addition, prehistoric extreme wave deposits that may be recovered from this coastal region (e.g., Claro et al., 2021) may also be associated with the Sulu megathrust rupture and tsunami.

It should be noted, however, that these eyewitness accounts pose uncertainties and are only available in certain coastal areas that have been inhabited before these events. For events before the 1960s, the estimated earthquake magnitudes are based on reported damages or anecdotes. The magnitude of the 1897 M_s 7.5, for example, may have been underestimated. This event may have also been a tsunami earthquake which generates large tsunami waves compared than its earthquake magnitude (e.g., Tanioka and Seno, 2001), or an earthquake-triggered submarine landslide also contributed to the basin-wide tsunamis as these submarine features are prevalent in this active margin (Nawanao and Ramos, 2023). As our models assumed blind megathrust rupture, these historical events may also have been triggered by other rupture mechanisms (e.g., outer rise, tsunami earthquake, intraplate, upper-plate earthquakes) with a different set of parameter estimations. With the scarcity of information, these historical accounts are still valuable to compare with tsunami models. What is consistent is that the coastal positions and configuration relative to the trench dictate where the highest exposures are identified, i.e., controlled by wave directivity.

5.3 Limitations and Future Work

A major limitation of this study is the uniform slip distribution along the megathrust fault plane, which is based on the average slip from the scaling equation. Earthquake slip distributions are typically heterogeneous, and they control whether the initial seafloor displacement can generate tsunamis. Megathrust ruptures with larger slip distributions downdip but lower on the updip can produce lower tsunami waves (e.g., Delouis et al., 2010; Fujii et al., 2020), which can lead to overestimation from these rupture scenarios. Conversely, large slip updip, especially for trench-breaching and splay faulting rupture types, can produce larger seafloor displacement (e.g., Baba et al., 2006; Fujiwara et al., 2011) that may underestimate the resulting tsunami wave heights from this study. Splay faulting and trench-breaching megathrust scenarios are also not accounted for in this study, as these modes of rupture are rarely reported, but they must be considered in future tsunami hazard assessments in this subduction zone. While the geometry of the Negros–Sulu subduction zone is currently not well constrained, the generally thick sediment cover overlying the Negros and Sulu subduction zone (Fig. 1) may decrease the slab dip and can widen the seismogenic zone (Brizzi et al., 2020). With the bathymetry data used in this study and the lack of variable and asperity-based slip distribution when considering megathrust ruptures, we recommend that the maximum wave heights are treated as minimum and that the actual highest tsunami waves are more likely underestimated.

Future studies along the Negros–Sulu subduction zone should aim to constrain both the geometry of the subduction slab and coupling ratios where slip deficits can be derived as applied in other subduction zones (e.g., Chlieh et al., 2008; Métois et al., 2013, 2016; Michel et al., 2019; Moreno et al., 2010; Ozawa et al., 2012). The role of submarine landslides in tsunami generation can also be explored in future studies as these features are prevalent in this active margin (Nawanao and Ramos, 2023). Submarine landslides have been shown to amplify tsunamis when the failure is coeval to megathrust earthquakes or can also be the main driver of localized but large tsunamis (e.g., Heinrich et al., 2001; Imamura et al., 1995; Nakata et al., 2020; Ramirez et al., 2022; Schambach et al., 2020; Takagi et al., 2019). Other mechanisms such as tsunami earthquakes and outer rise faulting may also be possible where focal mechanism solutions indicate normal faulting along these sections of the subduction zone (Fig. 2).

Updated seismic reflection and geophysical surveys across the Negros–Sulu Trench may also reveal underlying complexities that can influence the potential to generate megathrust rupture, including splay faults, fluids, and subducting seamounts. Gathering high-resolution bathymetry data is necessary to capture the local-scale wave amplification and deformation. The accuracy of nearshore bathymetry connecting inland and offshore must also be assessed, as the nearshore width shown in this study has significant effects on the resulting maximum tsunami waves along the coasts. Furthermore, mapping of pre-existing natural and man-made barriers for site-specific simulation and vulnerability assessment will also be important to assess their effectiveness in minimizing tsunami inundation.

The identified localities with the highest exposure to tsunami waves are potential sites (e.g., V, VI, XII, XIII) in searching for paleotsunami deposits that may serve as records for prehistoric earthquakes. Most importantly, the estimated moment magnitude, maximum tsunami wave heights, flow velocity, and minimum arrival times can provide invaluable information in creating strategies and policies for local government offices and coastal communities to mitigate hazard risk. These include information dissemination in coastal communities, delineation of evacuation routes and sites, and proper engineering designs such as coastal forests, green embankments, seawalls, offshore breakwaters, and dikes (e.g., Benazir et al., 2024; Dengler, 2005; Oetjen et al., 2022; Strusińska-Correia, 2017). In addition, as the majority of these coastal areas tend to have short arrival times, site-specific inundation modeling and geotechnical analysis are recommended for identifying and creating infrastructure for vertical evacuation.

6 Conclusion

We modeled the fault rupture and tsunami hazard potential of the Negros–Sulu subduction zone by constraining six megathrust fault segments at three dip angle scenarios (10, 20, and 30°). Rupture-scaling equations using the fault plane area estimate potential magnitude ranging from M_w 8.0 to 8.6 for the Negros Trench and M_w 8.1 to 8.9 for the Sulu Trench segments. From the estimated source parameters, elastic dislocation, and 4 h tsunami simulation, coastal areas directly facing the trench segments (VI Sipalay–Bayawan, XII Sindangan–Labason, XIII Siocon–Sibuco, XIV Basilan–Sulu, XV Moro Gulf) segments are potentially exposed to 5–8 m wave height and 5–7 m s⁻¹ flow velocity with a minimum arrival time of <2 min in majority of these regions, demonstrating the major control of wave directivity. Coastal areas parallel and on the rear side of the trench (X Palawan, XI Cagayancillo) are prone to tsunami wave height and flow velocity of up to 2–4 m and 2–4 m s⁻¹, respectively. Due to the cumulative effects of wave dissipation, interference, diffraction, and refraction, interior and rearmost islands are also exposed to <2 m tsunami wave height and <2 m s⁻¹ flow velocity. We also identified regions in the Sulu Sea where resonance on embayments can amplify tsunami waves, while nearshore width >30 km can lead to attenuation before reaching coastlines.

The multitude of future work reflects the scarcity of information along this active margin and the opportunity to fill these gaps. What this study provides, amid the limited information, are the maps of the minimum arrival, maximum tsunami wave, and flow velocity across the entire central and southern Philippines from different scenarios, informed based on available data and current understanding of megathrust ruptures. Amid the limitations of these models, we recommend that the tsunami wave heights and flow velocity for megathrust rupture scenarios are treated as minimum, serving as baseline reference in these coastal areas. This new information is useful not only for future studies to extend our records by searching paleotsunami deposits but most especially in conducting regional-scale mitigation assessments of coastal areas to the imminent threat of tsunamis that have severely affected the Philippines in the past.