Revised earthquake sources along Manila trench for tsunami hazard assessment in the South China Sea

Seismogenic tsunami hazard assessments are highly dependent on the reliability of earthquake source models. Here in a study of the Manila subduction zone (MSZ) system, we combine the geological characteristics of the subducting plate, geometry, and coupling state of the subduction interface to propose a series of fault rupture scenarios. We divide the subduction zone into three rupture segments: 14–16, 16–19, and 19–21.7^∘ N inferred from geological structures associated with the down-going Sunda plate. Each of these segments is capable of generating earthquakes of a magnitude between $M_{\mathrm{w}}=\mathrm{8.5}+$ and $M_{\mathrm{w}}=\mathrm{9}+$, assuming a 1000-year seismic return period as suggested by previous studies. The most poorly constrained segment of the MSZ lies between 19 and 21.7^∘ N, and here we use both local geological structures and characteristics of other subduction zone earthquakes around the world, to investigate the potential rupture characteristics of this segment. We consider multiple rupture modes for tsunamigenic earthquake and megathrust-splay fault earthquakes. These rupture models facilitate an improved understanding of the potential tsunami hazard in the South China Sea (SCS). Hydrodynamic simulations demonstrate that coastlines surrounding the SCS could be devastated by tsunami waves up to 10 m if large megathrust earthquakes occur in these segments. The regions most prone to these hazards include west Luzon of Philippines, southern Taiwan, southeastern China, central Vietnam, and Palawan Island.

1 Introduction

Large and damaging tsunamis are commonly triggered by megathrust ruptures that occur along convergent plate boundaries (i.e., subduction zones). Since 1900, many megathrust ruptures have triggered numerous devastating near- and far-field tsunamis including the 1952 M_w=8.8–9.0 Kamchatka event (e.g., Johnson and Satake, 1999; Kanamori, 1976), the 1960 M_w=9.5 event in the Chile subduction zones (e.g., Cifuentes, 1989; Moreno et al., 2009), the 1964 M_w=9.2 Alaska earthquake (e.g., Plafker, 1965), the 2004 M_w=9.2 Sumatra–Andaman earthquake along the northern Sunda Trench (e.g., Vigny et al., 2005; Banerjee et al., 2007; Chlieh et al., 2007), and the more recent 2010 M_w=8.8 Maule event in Chile (e.g., Vigny et al., 2011; Pollitz et al., 2011) and 2011 M_w=9.0 Tohoku-Oki earthquake along the northwest border of the Pacific Ocean (e.g., Koketsu et al., 2011; Wei et al., 2012). These earthquakes and their associated subduction zones have been intensively studied from different perspectives, including their tectonic settings and long-term evolution, seismic activities, and geodetic and geophysical features. In contrast, the Manila subduction zone (MSZ), which extends from southern Taiwan to the southern tip of Luzon Island in the Philippines along the eastern margin of the South China Sea (SCS) (Fig. 1), receives less attention, even though it shares many similarities with megathrust systems where large tsunamigenic earthquakes have occurred (Hsu et al., 2012, 2016).

Figure 1 Tectonic setting and historical tsunami catalogs in the South China Sea region. Colored circles indicate published tsunami catalogs and are labeled in the legend. Red triangles represent historical tsunamis related to volcanic activities. Red barbed curves show the megathrusts in this region. Geological tsunami records are marked with red stars (Ramos et al., 2017; Sun et al., 2013; Yang et al., 2018; Yu et al., 2009). The megacities are labeled in the legend and the seafloor subducting features are highlighted in the map. The historical earthquakes with M_w>6.5 in the Philippines are labeled. The likely tsunami events reported in the megacities are also labeled in the map. The two dashed lines represent the possible trace of the subducted Scarborough seamounts underneath the overriding plate as imaged from the tomography study (Wu et al., 2016). The rupture zones are denoted by the colored–shaded curves.

Over the past decade, attempts to study megathrust earthquakes and tsunamis from the Manila subduction zone have started to gain momentum. A number of rupture models have been used to assess potential tsunami hazard in the SCS (e.g., Hong Nguyen et al., 2014; Liu et al., 2009; Megawati et al., 2009; Okal et al., 2011; Wu and Huang, 2009), and yet the simulated tsunami wave heights and the subsequent hazard assessments differ greatly among studies (Hong Nguyen et al., 2014; Liu et al., 2009; Megawati et al., 2009; Okal et al., 2011; Wu and Huang, 2009; Xie et al., 2019). The difference often lies in the proposed fault-slip magnitudes of these models, and also the fault geometries used. Large variability in the results produced by these models underscores the fact that the seismogenic behaviors of the MSZ are still poorly understood. Some of the challenges which stand out and need to be resolved include assessing whether the MSZ is capable of hosting $M_{\mathrm{w}}=\mathrm{9}+$ earthquakes, investigating the amount of tectonic strain it has accumulated and its style of strain accumulation, and constraining how that strain is likely to be released in future.

Several lines of evidence suggest that the Manila trench has the potential to host a giant rupture capable of generating a basin-wide tsunami. Firstly, both historical earthquake records and modern seismicity databases (Hsu et al., 2012, 2016) indicate an absence of earthquakes larger than M_w=7.6 since Spanish colonization of Luzon in the 1560s (Bautista et al., 2012; Megawati et al., 2009; Ramos et al., 2017; Terry et al., 2017). The lack of significant megathrust-related earthquakes in modern records implies either a predominately aseismic megathrust or a highly coupled interseismic megathrust with the potential for a large ($M_{\mathrm{w}}=\mathrm{8.5}+$) rupture (e.g., Hsu et al., 2012). Several recent studies favor the high interseismic coupling model since both the analysis of earthquake focal mechanisms and geodetic monitoring results demonstrate that the upper plate is under shortening, which suggests that the megathrust, at least since the 1960s, shows minimal creeping (Bautista et al., 2001; Hsu et al., 2012, 2016). Secondly, the rate of plate convergence across the Manila trench is up to 90–100 mm yr⁻¹ – faster than the convergence rate of the Sumatra, Japan, and Nankai subduction zones, all of which have hosted giant earthquakes in the past few decades (McCaffrey, 2008; Megawati et al., 2009; Hsu et al., 2016, 2012;). Since the MSZ did not produce any significant events in the past 4 centuries, >30m of slip deficit is estimated to have been accumulated on the subducting interface (Megawati et al., 2009; Hsu et al., 2016). Thirdly, historical documents together with a few geological records across the SCS basin have reported nearly 130 tsunami events with different generation mechanisms (i.e., earthquakes, submarine landslides, volcanic eruptions). Although the credibility levels of these records varies (Bautista et al., 2012; Lau et al., 2010; Paris et al., 2014) and the geologically based interpretation suffers from the challenges of distinguishing tsunami waves from extreme storm surges, a series of records stand out with a similar range of event ages. Notably, four independent geological and geomorphological studies (Ramos et al., 2017; Sun et al., 2013; Yang et al., 2018; Yu et al., 2009) (Fig. 1) have purported evidence from coastal deposits, which they have inferred to be the result of a large tsunami event in SCS around 1000 to 1064 CE, which is nearly coincident with a large historical wave event recorded in Chao'an, Guangdong, in November, 1076 CE (Lau et al., 2010). The four independent sites of geological evidence are located at Dongdao island (Sun et al., 2013), Yongshu island (Yu et al., 2009), Badoc island near Luzon (Ramos et al., 2017), and Nan'ao Dao on the southern Chinese coastline (Yang et al., 2018) (Fig. 1). Since these studies identified only one event and if it were indeed generated by one tsunami, then we can conclude that the event was likely to be basin-wide and triggered by a very large MSZ event. Such an event would be a megathrust earthquake with sufficiently large rupture up to 1000 km long. Such a long and persistent rupture is comparable to the rupture length of the 2004 Sumatra–Andaman earthquake (e.g., Megawati et al., 2009). With the aforementioned evidence, there is no reason to rule out the possibility that the Manila trench could rupture as an M_w=9 earthquake (i.e., Megawati et al., 2009; Hsu et al., 2016). The current status of the Manila subduction zone could be an analog of the Sumatran subduction zone before the 2004 M_w=9.2 Sumatra–Andaman event between Myanmar and Aceh where a paucity of earthquakes > M_w=8 precede the 2004 event (Chlieh et al., 2008; Hsu et al., 2012), despite the very different geological settings (i.e., age, buoyancy, fault geometry) between these two subduction zones.

The SCS region is vulnerable to potential tsunami hazard. It covers an area of ca. 3.5 million km² (Terry et al., 2017), and is encircled by the coastlines of southeastern China, southern Taiwan, the western Philippines, eastern Vietnam, northern Borneo, and eastern Malaysia, forming a semi-enclosed basin (Fig. 1). The SCS coastline is one of the world's most densely populated with more than 80 million people living in the surrounding coastal cities (Terry et al., 2017). Many of these coastal cities serve as the economic centers and play pivotal roles in their respective countries' economic development. The coastline also hosts a very high density of major infrastructure (i.e., nuclear power plants, ports, airports). Data from the World Nuclear Association show that more than 10 nuclear power plants are currently in operation or about to start construction on the SCS coastline (http://www.world-nuclear.org/information-library, last access: 9 April 2019). Thus, if a large megathrust earthquake (e.g., M_w>9) were to occur within the SCS basin (Li et al., 2018), the impact would be amplified and much more devastating as the SCS is only about 1∕20 the size of the Indian Ocean. It is therefore crucial to provide physically based earthquake rupture models for a more realistic tsunami hazard assessment in the SCS region.

This study differs from previous studies (e.g., Hong Nguyen et al., 2014; Liu et al., 2009; Megawati et al., 2009; Okal et al., 2011; Wu and Huang, 2009) because we utilize a geodetic coupling model constrained by 17 years of GPS velocity measurements (Hsu et al., 2016) to propose a suite of better constrained physically based earthquake rupture scenarios. We also consider rupture segmentations constrained by the geological characteristics and the relief of the subducting Sunda plate. Scenario-based rupture models are different with the probabilistically based tsunami hazard assessments within which hundreds and thousands are implemented for rupture uncertainty estimates. Therefore, the probabilistic approaches (e.g., Li et al., 2016; Grezio et al., 2017) are often more complex to understand and implement than the scenario-based approaches. Here the proposed rupture models afford a physically based understanding of the tsunami hazard in the SCS. As a demonstration, we implement the rupture models to conduct hydrodynamic simulations to assess the tsunami characteristics along the coastlines of the SCS.

2 Refined possible rupture scenarios

Forecasting the extent and the slip distribution of earthquake ruptures is a challenging task. Before the 2004 M_w=9.2 Sumatra–Andaman earthquake (Chlieh et al., 2008), a M_w=9 earthquake had never been anticipated along the Sunda Trench, due to its oblique convergence orientation and seismically inactive feature (Satake and Atwater, 2007). Globally, the eventual ruptures of some unexpected fault locations keep surprising scientists (Bilek and Lay, 2018). We have seen partial ruptures of fully locked megathrusts (Konca et al., 2008; Qiu et al., 2016; Ruiz et al., 2014; Schurr et al., 2014), and piecemeal breaks in the center of perceived seismic gaps (e.g., Salman et al., 2017). Even with improved observations, it remains difficult to constrain the magnitude of a potential earthquake in the first order, and even more difficult to define the rupture pattern (e.g., Lay, 2018). A recent example comes from Japan where Loveless and Meade (2010) used a number of inland GPS stations to estimate the coupling state of the Japan megathrust before the 2011 Tohoku-Oki earthquake. They indicated the spatial extent of a possible future rupture. Notably, the rupture models constrained by multiple geodetic data sets after the 2011 earthquake (Koketsu et al., 2011; Loveless and Meade, 2011; Wei et al., 2012) are significantly different to the coupling map of Loveless and Meade (2010). The discrepancy between a coupling map and actual rupture estimates has also been observed at other subduction zones (e.g., Ruiz et al., 2014; Schurr et al., 2014) and for the collision zone between the Indian and Eurasian plates (Avouac et al., 2015; Qiu et al., 2016; Stevens and Avouac, 2015). Clearly, our current knowledge of the seismogenic characteristics of giant earthquakes remains deficient.

Great efforts have been made to investigate the physical parameters that characterize subduction zones with regard to the geometry, geology, and dynamics (Schellart and Rawlinson, 2013). Systematic analysis of collections of great earthquakes globally indeed suggests that some of the physical parameters do play key roles in controlling the rupture characteristics (Bilek and Lay, 2018; Bletery et al., 2016; Schellart and Rawlinson, 2013), although limitations in the historical earthquake records inevitably make it difficult to have high confidence on such relationships. Taking into account the geometrical effects, previous studies have divided the entire Manila subduction zone into three segments (i.e., Zhu et al., 2013; Li et al., 2016; Gao et al., 2018). Here we follow the segments proposed by Li et al. (2016), and we provide new constraints on earthquake and tsunami potentials by combining geological information and the geodetic constrained coupling map to adjust these segments accordingly. The modulated three segments are 14–16, 16–19, and 19–22^∘ N. Their significance is detailed in subsequent sections.

2.1 Rupture segment 1 (zone 1, 14–16^∘ N)

The Manila trench primarily starts from ca. 13^∘ N west of Mindoro and ends at ca. 22^∘ N southwest of Taiwan, and beyond these bounds the Manila trench gradually transforms into a collision and accretionary belt in the north and south (Fig. 1). At the southernmost area of the Manila trench, the strike direction of the trench bends to the southeast offshore of Mindoro Island (ca. 13^∘ N) before it further collides with Panay (ca. 11^∘ N). Within this region (ca. 13 to 11^∘ N) the relocated seismicity suggests the subducting slab dips almost vertically, with an absence of the deep seismicity (Bautista et al., 2001). Based on these features, Bautista et al. (2001) suggest the subducting slab may have been heated up and assimilated into the mantle. We, therefore, interpret that the great megathrust earthquake is less likely south of 13^∘ N. Li et al. (2016) placed the southern boundary of the first segment at ca. 12.5^∘ N. In contrast, Bautista et al. (2001) proposed a slab tear at ca. 14^∘ N, which is the result of the collision of a micro-continental plate with Mindoro and Panay islands and as evidenced by the narrow seismicity gap north of 14^∘ N that trends northeastward (Fig. 2). Based on these geological characteristics and geodetic measurements, together with the fact that the spatial coverage of GPS measurements in this region only allows us to estimate the coupling status starting at 14^∘ N to the north (Hsu et al. 2016), we move the southern boundary of the first segment from ca. 12.5^∘ N proposed by Li et al. (2016) to 14^∘ N, but we do not rule out the possibility of ruptures that propagate across 14 to 13^∘ N or even beyond.

Figure 2 Seismicity (M_w>4.5) in the Manila subduction zone between 1900 and 2018. This data set is downloaded from the USGS catalog. Color represents the depth and size scales of the seismic moment magnitude indicated in the legend.

Moving to the north, between 16 and ca. 17.5^∘ N, a bathymetric high called the Scarborough seamount chain is subducting beneath the Philippine plate. The Scarborough seamount can be traced between ca. 12 and 18^∘ N from the subducting Sunda plate and between ca. 16 and 19^∘ N (Fig. 1) after subducting beneath the Philippine plate from the regional tomography model (Wu et al., 2016). This seamount chain has been interpreted as part of an extinct middle ocean ridge (MOR) that is either presently being accreted or subducted under the trench at 16 ^∘N (Ludwig, 1970; Pautot and Rangin, 1989). A slab tear was proposed at 16^∘ N based on seismic-related strain energy release of intermediate depth and shape changes in the dip angle of the slab (Bautista et al., 2001). Although great earthquakes can rupture across the seamounts or morphological bounds occasionally (e.g., Bell et al., 2014; Duan, 2012; Kumagai et al., 2012), global observations suggest that in many cases seamounts or barriers impede (Singh et al., 2011; Wang and Bilek, 2011) or confine rupture propagations (Qiu et al., 2016). Further, we note that slab tears at 14 and 16^∘ N bound the southern and northern tips of the highly coupled west Luzon trough (Hsu et al., 2012, 2016) coincidently, and these tears may act as morphological barriers to limit the rupture propagation similar to that noted from the 2015 M_w=7.8 Nepal event (Qiu et al., 2016). We, therefore, define the region between 14 and 16^∘ N as segment 1 (zone 1) (Fig. 3a and d).

Figure 3 Proposed rupture slip models based on coupling models from Hsu et al. (2016) assuming a 1000-year seismic return time period. Panels (a) and (b) show the slip models from coupling model A (Hsu et al., 2016) in zones 1 and 2, respectively. Panel (c) shows a proposed hybrid model based on coupling model A (19 to 20^∘ N) and a Gaussian shape of slip distribution (20 to 21.7^∘ N) with a 50 % coupling ratio in zone 3. Panels (d–f) represent the same slip models as (a–c) but based on coupling model B (Hsu et al., 2016). CM refers to coupling model. Coupling models A and B are from Hsu et al. (2016) that are shown in the inset map. White arrows show the possible slip directions during an earthquake. Vectors in the coupling maps show the slip deficit direction that is accumulated for future release in earthquakes. The estimated seismic moments of each model are labeled in each subplot with rigidity of 30 GPa. The slip magnitude and coupling ratio are shown by its corresponding color scales.

2.2 Rupture segment 2 (zone 2, 16–19^∘ N)

As noted in Sect. 2.1, the Scarborough seamount chain is located between ca. 16 and 17.5^∘ N where the subducting Sunda plate meets with the Philippine plate (Fig. 1). A regional tomography model also suggests that the subducted seamount chain can be traced between ca. 16 and 19^∘ N (Wu et al., 2016). In this subducted seamount region, the absence of seismicity and seismic-related strain energy release at intermediate depths suggest the possible trajectory of the MOR that is interpreted to be still hot and deforming plastically (Bautista et al., 2001). Globally studies of subducting seamount systems suggest that large fracture zones are formed surrounding the seamount, and the highly fractured region can act as barriers to hinder the rupture propagation (e.g., Wang and Bilek, 2011). Because the stress concentration in and around the fracture zones is high and may easily reach failure criteria, the seamount can trigger (e.g., Kumagai et al., 2012; Koyama et al., 2013) the failure of highly stressed asperities in the neighborhood, nucleating as a great earthquake (e.g., Kumagai et al., 2012; Koyama et al., 2013). Previous studies also suggest that seamounts cause persistent fault creep (e.g., Singh et al., 2011) or rupture as small earthquakes due to localized areas of high fracture and associated regional stress anomalies (e.g., Wang and Bilek, 2011). Thus, fault creep and the rupture of single or multiple asperities are all possible in this region.

The Geodetic coupling map constrained by long-term GPS velocity measurements indicates that the seamount chain region (i.e., ca. 16 to 19^∘ N) is less coupled (Fig. 3, coupling models A and B), partially due to the fault creep caused by the seamounts or poor constraints by paucity of the offshore observations (Hsu et al., 2012, 2016). The weak coupling extends further north to 19^∘ N, in the area of the southern tip of the North Luzon Trough and west of the northern tip of Luzon Island. This area is likely creeping or weakly coupled (Fig. 3, coupling mode A and B). Additionally a trench-parallel gravity anomaly (TPGA) has been interpreted with great subduction earthquakes occurring predominately in areas characterized by strongly negative TPGA, while regions with strongly positive TPGA are relatively aseismic (Song and Simons, 2003). We note that positive TPGA covers from ca. 16 to 19^∘ N (Hsu et al., 2012), coinciding with the geodetically determined weakly coupled and creep regions. Considering all these factors mentioned above, we redefine segment 2 (zone 2) as the region between 16 and 19^∘ N as (Fig. 3b and e) slightly extends further north when compared with the same segment of Li et al. (2016).

2.3 Rupture segment 3 (zone 3, 19–22^∘ N)

The area of the megathrust bounded between the southern tip of Taiwan and northern Luzon (between 19 and 22^∘ N) (Fig. 1) is poorly understood, as the current available geodetic measurements are sparse and primarily deployed in the volcanic islands to the east which are far away from the Manila trench (Hsu et al., 2012, 2016). In this region, the Manila trench bends sharply at 20^∘ N (Fig. 1). Geologically the bending has been interpreted as the result of the subduction of a high-relief bathymetrical plateau that is sufficient buoyant to impede subduction (Bautista et al., 2001; Suppe, 1988) or may due to thick sediments (Lin et al., 2009). Additionally, here regional block faulting stretches the continental crust, resulting in numerous micro-continental fragments. Further, the 1980s geophysical studies (Taylor and Hayes, 2013) have recovered a magnetic quiet zone characterized to the continental-to-oceanic boundary (Bautista et al., 2001), and this zone was further interpreted with a transition zone between a continental and oceanic lithosphere (Taylor and Hayes, 2013). If these numerous fragments are indeed subducting beneath the Philippine Sea Plate, then they would have to be buoyant enough to resist the subducting process at 20^∘ N with fast subducting of the neighboring portions of the trench that may extending south to 19^∘ N. Such a situation would result a complex stress field in the upper plates that were mirrored by diverse and complicated focal mechanism solutions (Bautista et al., 2001).

As more marine geophysical data becomes available, there is an increased understanding of the geological structure and potential seismogenic faults (Lin et al., 2009). Detailed analysis of seismic reflection data (i.e., line 973 in Lin et al., 2009) reveals prominent seismogenic structures in the region, which include frontal decollement beneath the lower-slip domain and out-of-sequence thrusts (OOSTs) in lower- and upper-slope domains (Lin et al., 2009; Zhu et al., 2013). Evidence from the thermal regime of these structures suggests that the megathrust and part of the frontal decollement are seismogenic (Lin et al., 2009). These seismogenic structures are found to be analogous to that observed in the Nankai prism of the Nankai Trough, Japan, posing potential for generating great earthquakes and tsunamis as they did in Nankai (Lin et al., 2009; Yokota et al., 2016).

Fan et al. (2016) revealed a low-velocity zone that spans from shallow to deep depths of 20–200 km beneath the prism, suggesting that the collision develops northward and the subducting process may stop at 22^∘ N. Coincidently, at a similar latitude (21.5^∘ N), Lin et al. (2009) interpret that south of 21.5^∘ N, the subducting is active while north of this latitude the plate convergence is accommodated by intense compressional deformation of the crust due to the buoyance of the Eurasian plate that resists subduction. Consequently, in light of the geological evidence noted above, we slightly shorten the northern boundary of segment 3 from Li et al. (2016), and we define the region between 19 and 21.7^∘ N as segment 3 (zone 3) (Fig. 3c and f).

3 Proposed slip deficit models

Using geodetic surface measurements, a velocity value can be derived and used to constrain the elastic strain accumulation rate between the subduction plate interfaces, the so-called interseismic coupling model (Chlieh et al., 2008; Hsu et al., 2012, 2016; Loveless and Meade, 2010; Megawati et al., 2009). This model reveals strain accumulation within seismic cycles that can potentially be released during great earthquakes, although the final rupture extent is commonly not exactly the same as forecasted by the coupling maps (Konca et al., 2008; Ruiz et al., 2014) and in some cases uncoupled parts of the megathrust may regularly produce tsunamis (Witter et al., 2016). However to move towards an associated tsunami hazard assessment from such potential ruptures, the coupling map, although not perfect, is often the necessary choice (e.g., Power et al., 2012; Megawati et al., 2009).

Using decades-long GPS velocity measurements, Hsu et al. (2016) proposed two coupling models (A and B) that best explain the plate movements and coupling state on the Manila megathrust and other faults on Luzon Island. With this coupling or slip deficit rate estimates and the possible seismic return time period, we can forecast the likely slip distributions that may fail in future earthquakes. For seismic return time period, given the short duration of historical records relative to the return periods of large-magnitude events of interest, and limitations in our capacity to infer earthquake return periods from first-principle physics, it is unrealistic to expect to develop a comprehensive understanding of seismic return periods. We thus have to rely on the observations. The modern seismic records for the Manila trench only trace back to ∼1900 and provide constraints on the natural frequency of earthquakes with its corresponding magnitude assuming the Gutenberg–Richter (G–R) earthquake relations, and thus are often implemented for tsunami hazard assessment (Li et al., 2016; Power et al., 2012). Historical records since the 1560s suggest that there is no recorded earthquake with M_w>7.6 in the Manila subduction zone, implying that the determined return time period for a great earthquake from the G–R relation will likely be poorly constrained (Hsu et al., 2016). However geological evidence from purported tsunami deposits may provide evidence of tsunamis at four locations in the SCS (i.e., Fig. 1, Ramos et al., 2017; Sun et al., 2013; Yu et al., 2009; Yang et al., 2018). Some studies suggest that a giant tsunami event might have occurred ca. 1000–1064 CE (Ramos et al., 2017; Yang et al., 2018). With an assumption of a 1000-year return period, the magnitude can reach $M_{\mathrm{w}}=\mathrm{9}+$ from geodetic analysis (Hsu et al., 2016). Here we choose to model scenarios that release 1000 years of accumulated strain because these represent large, rare, and yet plausible events, which are of interest for hazard assessment purposes, and paleogeologic data indicate that large events may have occurred about 1000 years ago.

Based on coupling models A and B of Hsu et al. (2016) in which the spatial distribution of slip rate and coupling rate is available, we use a return period of 1000 years to calculate the slip deficit of great earthquakes assuming each event releases 1000 years of strain accumulation while ignoring the possible portion of strain released by smaller events. For the predefined zones 1 to 3 (see Sect. 2.1–2.3), different approaches are used. For zones 1 and 2 where the coupling ratios and slip rate are relatively better constrained than zone 3, we calculate the slip deficit by multiplying the slip deficit rate at each triangle node (Fig. 3a, b, d, and e) with 1000 years. The slip deficit models in zone 1 for models A and B (Fig. 3a and d) are similar with the maximum slip > 50 m occurring at ca. 20–30 km seismogenic depth due to the high coupling ratio. For zone 2, the slip model based on A has a compact area and less slip compared with the slip model based on B (Fig. 3b and e). This is because the extra North Luzon Trough fault was introduced in model B, resulting in larger spatial extent and higher coupling while equally explaining the GPS velocity measurements (Hsu et al., 2016).

Due to paucity of observations in zone 3, no coupling ratios were resolved. Geologically this zone is much more complicated than zones 1 and 2 (Lin et al., 2009). Multiple OOSTs are revealed from seismic reflection profiles (Lin et al., 2009; Zhu et al., 2013). Failure of these OOST (also called megasplay) faults with high dip angle contributes to generating devastating waves as evidenced by historic tsunami events in other subduction zones (Moore et al., 2007; Park et al., 2002). It is, therefore, crucial but difficult to precisely quantify the individual role of the OOSTs and megathrust in tsunami generation.

We propose two end-member scenarios, considering different rupture modes in zone 3 with two steps. We first calculate the slip deficit from the slip deficit rate of models A and B between 19 and 20^∘ N. We then consider two end-member scenarios in the region from 2 to 21.7^∘ N. The first member is the seismogenic events with rupture depths determined from a collection of global centroid moment tensor (GCMT) solutions of the world megathrust earthquakes (Fig. 4). We assume the fault slip pattern follows a Gaussian distribution centered at 25 km of the mean depth from the global great earthquakes. We cut off slip deeper than 50 km as the rock properties at this depth and beyond induce semi-brittle and ductile flow (Hippchen and Hyndman, 2008; Hyndman and Wang, 1993; Wang, 2007). This can capture to first order the potential slip extent (Fig. 3c and f), with a depth range of slip consistent with observations from global megathrust great earthquakes (e.g., Chlieh et al., 2007; Pollitz et al., 2011; Ruiz et al., 2014; Salman et al., 2017; Wei et al., 2012). For the second mode, we consider tsunamigenic events similar to the 2011 M_w=9.0 Japan earthquake in which the earthquake can rupture all the way to the trench. We estimate the plate convergence rate in the fore-arc in zone 3 to be 67 mm yr⁻¹ (Hsu et al., 2009) with a 24.5 mm yr⁻¹ shortening under the 91.5 mm yr⁻¹ plate convergence rate with respect to the Sunda plate (Hsu et al., 2016; Sella et al., 2002). We assume 67 mm yr⁻¹ convergence was fully accommodated by the megathrust and implement it as the amplitude of the Gaussian distribution, allowing the maximum slip occurring at the trench (Fig. 5a and b). For each rupture mode, we have two slip models corresponding to coupling model A and B, and assume half of plate convergence rate is accommodated by the megathrust (Fig. 5a and b, with 80 % coupling ratio shown in Fig. 6c and d). For the second-member model, we implement rupture on both the megasplay fault and the megathrust assuming each of them to accommodate half of the fore-arc plate convergence and a uniform slip on the splay fault as a simple case (Fig. 5c and d). We implement this splay fault only with seismogenic rupture events as we think this case is easier due to the splay fault's bottom cut to the megathrust at seismogenic depth (Lin et al., 2009). We consider a 50 % coupling ratio for both the megathrust and splay fault (Fig. 5c and d, with 80 % coupling shown in Fig. 6a and b). Details about these proposed rupture scenarios are given in the summary Table S1 in the Supplement.

Figure 4 Depth distribution of the seismicity in the Manila subduction zone between 1970 and 2018. This data set is downloaded from the GCMT catalog. Colors represent the seismic moment magnitude. The giant 2004 Sumatra and 2011 Japan earthquakes are highlighted in the map.

Figure 5 Proposed slip models in zone 3. Panels (a, b) show the shallow rupture type of slip models (e.g., Tohoku-Oki) based on coupling models A and B (Hsu et al., 2016), respectively. Panels (c, d) represent megathrust (Fig. 3c and f) rupture and the out-of-sequence megasplay type of slip models, respectively. We assume 50 % coupling for the megathrust and the megasplay faults. CM refers to the coupling model shown in the inset map. White arrows show the possible slip directions during earthquakes. Vectors in the coupling maps show the slip deficit direction that is accumulated for future release in earthquakes.

Figure 6 (a) Seismogenic megathrust rupture together with megasplay rupture scenario with 80 % coupling ratio for each of them from model A of Hsu et al. (2016). Panel (b) same as (a) but from model B of Hsu et al. (2016). (c) Shallow rupture (e.g., Tohoku-Oki rupture) the same as Fig. 3a but with 80 % coupling ratio on the megathrust. (d) Shallow rupture (e.g., Tohoku-Oki rupture) same as Fig. 3b but with 80 % coupling ratio on the megathrust. A 1000-year seismic return period was assumed in the slip calculation.

The geometry of the OOST is derived from Lin et al. (2009) and covers the area from 20 to ca. 22.2^∘ N, as we ignored the bending portions of the OOST in the north and south although they still can rupture with a low probability. The fault is ca. 260 km long and ca. 16 km wide and it strikes 345^∘ to the north and dips 50^∘ to the east (Figs. 1, 5, and 6).

4 Tsunami impacts in SCS

4.1 Tsunami simulation setup

We use the Cornel Multi-grid Coupled Tsunami (COMCOT) model to simulate the hydrodynamic process of the tsunami waves (e.g., Wang et al., 2008; Philip, 1994; Li et al., 2016, 2018) produced by those proposed earthquake ruptures. The initial surface elevations generated by all the proposed rupture models can be found in the Supplement. To account for the nonlinear effect in the nearshore region, the simulation solves the nonlinear shallow water equations in spherical coordinates for the entire SCS region with a bottom Manning friction coefficient of 0.013 (Li et al., 2018). We used the 1 arcmin grid of General Bathymetric Chart of the Oceans (GEBCO) data for the modeling. A uniform grid was used because we do not focus on near- and onshore processes where high-resolution topographical data and good understanding of the bottom friction effect are required. Synthetic gauges along the 20 m isobaths are specified to record the tsunami waveforms. For the initial tsunami waves, we assume the rupture occurs instantaneously and the vertical seafloor deformation produced by the ruptures is equal to the initial ocean surface deformation (e.g., Li et al., 2016, 2018; Liu et al., 2009).

4.2 Maximum tsunami wave height

For all the simulated scenarios, the resulting wave height in the near-source regions mainly depends on the rupture location and earthquake magnitude. While in the relatively far field, the tsunami wave directivity effects and bathymetry effects also play important roles (Figs. 7 and 8). We describe the tsunami impact of each pair of source models from south (zone 1) to north (zone 3). Slip models in zone 1 generate the largest tsunami waves (>10 m) in western Luzon (Fig. 7a and b). Central Vietnam experiences a similar tsunami height (4–8 m) to the intermediate far-field area, western Palawan. Southeastern China and southern Taiwan could be hit by up to 5 m tsunami waves (Fig. 7a and b). Moving to zone 2, the slip models show the significant difference in terms of both magnitude and slip distribution between models A and B (Fig. 3b and e). Consequently, the tsunami impact caused by model B is much larger than that caused by model A in both near-source (e.g., western Luzon and southern Taiwan) and far-field regions (e.g., southeastern China and central Vietnam). Compared with the region most affected by slip models in zone 1, the worst-hit region also moves northward with the rupture location. Similarly, when the earthquakes rupture the megathrust in zone 3, the hardest-hit regions move further to the northern part of the SCS and concentrate in northern Luzon, southern Taiwan, and southeastern China (Fig. 7e and f). Further, Fig. 8 shows the diverse tsunami impacts generated by rupture scenarios in zone 3. Not surprisingly, the results suggest rupture models with higher coupling cases (Fig. 8a and b) result in larger tsunami wave heights in regions located in the northeast SCS despite the tsunami generation efficacy of shallow slip earthquakes (Fig. 8c and d). One interesting phenomenon worthy of mention is the high tsunami hazard of southeastern China regardless of the rupture locations. This is likely explained by the combined effect of tsunami wave directivity and bathymetry (Figs. 7 and 8). Tsunami waves refract significantly on the southern Chinese coast due to the shape and gradient of the continental slope, leaving southeastern China (including coastlines of Guangdong, Hong Kong, and Macau) in the direct tsunami path.

Figure 7 Modeled maximum tsunami wave heights and arrival time contours in the SCS. Panels (a, c, e) show the maximum tsunami wave heights generated from rupture zones 1–3 based on slip models calculated from coupling model A (Hsu et al., 2016), respectively. Panels (b, d, f) show the same maximum tsunami wave heights but with slip models calculated from coupling model B (Hsu et al., 2016). In zone 3, we show Gaussian slip distribution with a 50 % coupling ratio scenario with other example scenarios shown in Fig. 8. The solid black contours show hourly tsunami arrival time with half-hour increments (dashed contours). The colored dots show the subsampled location at 20 m water depth, with color showing the maximum wave heights.

Figure 8 Maximum tsunami wave heights from different rupture characteristics in zone 3 with hybrid coupling models. Panels (a, c, e) show the maximum tsunami wave heights based on coupling model A (Hsu et al., 2016). Panels (b, d, f) show the same maximum tsunami wave heights but with slip models based on coupling model B (Hsu et al., 2016).

To summarize, the near-source regions including western Luzon, northern Luzon, and southern Taiwan face the greatest tsunami hazard. The second most threatened areas are southeastern China, central Vietnam, and western Palawan. Archipelagos inside the SCS including Dongsha, Zhongsha, and Xisha also suffer severe tsunamis (up to 6–8 m tsunami wave height) when large earthquakes occur in zones 2 and 3. Coastal regions of northern Borneo, eastern Malaysia, eastern Thailand, and southern Cambodia are significantly less affected.

4.3 Tsunami travel time

The tsunami travel time is key information in tsunami evacuation planning. Similar to the other subduction zones, the near-source areas including the coast of Luzon and southern Taiwan suffer the highest tsunami waves with least evacuation time (Figs. 7 and 8). We plot the time series of tsunami waves generated by all the source models in selected synthetic gauges near nine major coastal cities in Fig. 9. Depending on the rupture locations, the tsunami arrival time is in minutes or less than half an hour for near-source cities, like Vigan, Kenting, and Kaohsiung (Fig. 9), posing great challenges to the early-warning system and subsequent evacuation process. In other areas tsunami wave travel time is relatively longer for Vietnam and southeastern China. The arrival time is commonly between 2 and 3 h after the earthquake for central Vietnam and 3–4 h for southern China. For the archipelagos inside the SCS, the tsunami waves arrive much earlier than they do on the mainland in Vietnam and China, typically ∼1 h earlier. The earlier arrival time in archipelagos make them ideal locations for installing tsunami monitoring instruments (e.g., tide gauges or GPS; see Peng et al., 2019). Such measurements may provide timely constraints on wave height for the evacuations in far-field areas. Detailed inundation maps of the main coastal cities in this region are highly recommended for designing evacuation routes.

Figure 9 Simulated tsunami wave time histories at example coastal cities in the SCS region. Top, middle, and bottom of each subpanel show simulated waves from ruptures in zones 1–3, respectively. A and B represent coupling models A and B from Hsu et al. (2016). For rupture zone 3, we show the Gaussian slip distribution with 50 % coupling ratio cases (Fig. 3c and f) as examples.

5 Discussion

How and where earthquake rupture will occur on a plate boundary is challenging to forecast (Bilek and Lay, 2018; Satake and Atwater, 2007). A comprehensive understanding of a single megathrust behavior may be impractical since the seismic cycle is typically on the order of hundreds and thousands of years, much longer than instrumental records. Conversely, understanding megathrust behaviors over different subduction zones at different time stages of their cycle offers insights into rupture style and characteristics. Previous studies have intensively investigated giant subduction zone earthquakes, gaining useful insights into physical parameters that are related to developing giant ruptures. Such physical parameters include the subducting plate age, rate, and buoyance of the slab (Kanamori, 2006; Nishikawa and Ide, 2014; Ruff and Kanamori, 1980, 1983); the forearc structures (Song and Simons, 2003; Wells et al., 2003); upper plate characteristics including plate motion (Schellart and Rawlinson, 2013); trench characteristics of the long-term migration (Schellart and Rawlinson, 2013) and sediment thickness (Heuret et al., 2012); and the width of seismogenic zones (Hayes et al., 2012; Schellart and Rawlinson, 2013). As more data become available, we understand that the convergence rate does not play a significant role in generating giant earthquakes (McCaffrey, 1994; Stein and Okal, 2007; Nishikawa and Ide, 2014). The maximum moment magnitude of a potential earthquake is often determined from seismic catalogue data, alternatively determined from basic moment conservation principles and catalog data (Rong et al., 2014; Kagan and Jackson, 2013). Overall, with current short observation time span compared with multi-century seismic return period, it is improper to make a determination on the relationship between these physical parameters and how big or how often a giant earthquake can occur in any subduction zone (McCaffrey, 2008). Clearly, long-term and complete observations within seismic cycles are required for a better understanding of subduction zone rupture behaviors.

Recently a summary study based on global subduction zone observations concluded that megaseismic events preferentially rupture flat, gently dipping interfaces (Bletery et al., 2016). In the Manila trench, the dip is gentle and progressively increases from north to south (Bautista et al., 2001). In zone 3, the presence of subducting plateau of the continental fragments results in a gently dipping, near-flat interface that potentially favors the development of giant earthquakes (Figs. 1 and 2). The dipping degree is in a similar range as those found in other subduction zones, e.g., Japan–Kuril–Kamchatka, Alaska–Aleutians, Sumatra–Java, South America, and Cascadia, which are known to produce M_w>9.0 earthquakes (Bletery et al., 2016).

Morphological barriers have been found to have a predominant role in controlling rupture propagation and style. The barriers can confine and arrest rupture propagation (Qiu et al., 2016), and act as a persistent fence to stop rupture (Meltzner et al., 2012; Morgan et al., 2017). Fault bends can also hinder rupture overstep at bending points (Wesnousky, 1988, 2006). In the case of the Manila subduction zone, the presence of the Scarborough seamount chain in zone 2 and slab tear in zone 1 indicates that a giant rupture propagation through zones 3 to 1 is less likely, although we do acknowledge that a rupture-across-zone earthquake is possible with very low probability, like the 2007 M_w=8.1 earthquake that ruptured a triple junction (Furlong et al., 2009; Taylor et al., 2008). Dynamic simulations do show possible scenarios that involve multiple portions of the Manila trench rupturing as a single giant earthquake (Yu et al., 2018). However, the details of the slab tear in zone 1 and the seamount chain in zone 2 were smoothed out in the simulation, due to the challenges of the numerical calculation (Yu et al., 2018).

Regarding the potential source of the geological records, the tsunami simulations suggest the difficulty of creating a scenario which could affect all four tsunami deposit locations with sufficiently high tsunami waves, especially for the record located on Yongshu island (Yu et al., 2009). Assuming all four records are indeed tsunami deposits, the spatial distribution demands the whole trench to rupture at once and the southern segment needs to extend further to 13^∘ in order to generate tsunami wave propagation in the southwest direction towards Yongshu. Another alternative explanation could be that the deposits on Yongshu island were generated by a large storm event instead of a tsunami event.

In summary, our definition of the rupture zones 1–3 are derived by taking into account the bathymetry features of the subduction Eurasian plate, earthquake focal mechanism distributions, structure-controlled TPGA, and more than 20-year-long GPS measurements. The refined coupling models (Hsu et al., 2016) offer more detailed images that reflect the likely motions on the plate interface. Combination of the coupling models and morphological bounds' constrained zone definitions provides more realistic rupture scenarios than planar fault with assumed uniform slip rupture cases. We have seen that finite rupture models of historical earthquakes indicate that slip is heterogeneous, and this is represented by our scenarios. Further detailed tsunami hazard assessment in the SCS demonstrates that uniform slip models underpredict tsunami hazards compared to a heterogeneous slip model (Li et al., 2016). Therefore, our refined earthquake rupture scenarios in zones 1 and 2 provide new insights into tsunami hazard assessment in the SCS. For zone 3, the scarcity of measurements and the presence of complicated geological structures result in a poor understanding of the seismogenic characteristics, although the tsunami-genic potential remains high (Lin et al., 2009). The possible ruptures provided in this study can be a first-order approximation of the earthquake scenarios in the region. Subsequent measurements collected in coming years can help us to refine our understanding in this region.

6 Conclusion

We have proposed updated earthquake rupture scenarios along the Manila trench based on new geological earthquake focal mechanism information and geodetic observations. These rupture models enable tsunami assessment in the SCS and subsequent detailed examination of inundation processes for megacities along the coastlines of the SCS.

Tsunami simulations based on these rupture scenarios indicate that the coastlines of the SCS region are under the risk of devastating tsunami waves, specifically for western Luzon of the Philippines, southern Taiwan, southeastern China, central Vietnam, and Palawan Island. In addition to the near-source region, southeastern China will also be attacked severely due to the bathymetry focusing effect no matter which portion of the Manila thrust breaks. Southern Taiwan is affected by ruptures in zones 2 and 3, with west Luzon affected by all earthquake scenarios. Central Vietnam and Palawan Island are mostly affected by ruptures in zones 1 and 2. In all cases, the waves sweep these coastlines within ca. 3 h. Our results highlight that it is necessary to conduct further detailed inundation investigations at these severely affected coastal regions, for future preparation on hazard mitigation plans. Our findings also provide useful information that could be used to find possible archived geological recordings of historical tsunami deposits, and call for subsequent paleo-sedimentology studies in the SCS basin.