Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 943147 10.3389/feart.2022.943147 Earth Science Original Research Bathymetric Highs Control the Along-Strike Variations of the Manila Trench: 2D Numerical Modeling Ma et al. Bathymetric Highs Subduction Ma Letian 1 2 * Chen Lin 3 Cheng Zihua 4 Gerya Taras 5 Li Jiabiao 1 2 1 Key Laboratory of Submarine Geosciences and Second Institute of Oceanography, Ministry of Natural Resources, Hangzhou, China 2 Southern Marine Science and Engineering Guangdong Laboratory, Zhuhai, China 3 State Key Laboratory of Lithospheric Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing, China 4 Key Laboratory of Ocean and Marginal Sea Geology, South China Sea Institute of Oceanology, Chinese Academy of Sciences, Guangzhou, China 5 Department of Earth Sciences, ETH Zurich, Zurich, Switzerland

Edited by: Rafael Almeida, San Diego State University, United States

Reviewed by: Fucheng Li, South China Sea Institute of Oceanology (CAS), China

Liming Dai, OUC, China

 *Correspondence: Letian Ma, malt@sio.org.cn

This article was submitted to Structural Geology and Tectonics, a section of the journal Frontiers in Earth Science

 11 07 2022 2022 10 943147 13 05 2022 17 06 2022 Copyright © 2022 Ma, Chen, Cheng, Gerya and Li. 2022 Ma, Chen, Cheng, Gerya and Li

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.

The Manila Trench is located at the eastern boundary of the South China Sea (SCS). It develops through the subduction of the SCS beneath the Philippine Sea Plate (PSP) since the early Neogene, driven by the northwestern plate motion of the PSP. The northern segment of the Manila trench at around 18° N—21.5°N is characterized by an obvious eastward convex in the trench shape and abrupt changes of slab dip angle, whereas the southern segment of the Manila trench at around 15°N—18°N is featured by an almost straight NS-trending trench line and smooth subducting slab morphology. However, the cause for the along-strike variations along the Manila trench remains poorly understood. In this study, we use 2-D thermo-mechanical modeling to investigate how bathymetric highs embedded in the subducting slab affect the topography of overriding plate and the morphology of subducting plate. Three major factors of bathymetric highs are systematically examined: 1) the crustal properties, 2) the width, and 3) the thickness. Geodynamic results suggest that the most important factor controlling abrupt changes in dipping angle is the crustal properties of bathymetric highs. Also, reduction of crustal thickness and increasing the width of continental bathymetric highs favor the abrupt change of dipping angle, whereas thicker (≥25 km) bathymetric highs are more likely to be blocked in the subduction zone before slab break-off. According to our numerical modeling results, we suggest that dramatic changes in the dip angle in the northern Manila trench and the convex shape were caused by subduction of a large thin continental terrane, whereas the smooth morphology of subducting slab in the southern segment and straight trench were associated with normal oceanic subduction with small seamounts.

 numerical modeling bathymetric highs subduction Manila trench South China Sea Introduction

Sharp discontinuities in trench shape and various morphologies of subducting slabs are often linked with bathymetric highs subduction (Miller et al., 2004, 2005; Mason et al., 2010). For example, the subduction of Ogasawara Plateau may have influenced the varied trench and slab morphology at the junction of the Izu-Bonin and Mariana arcs in the West Pacific (Mason et al., 2010). Similar arcuate plate boundary development is also found during the subduction of the Nazca Ridge beneath the continental South American Plate (Rosenbaum et al., 2005). These bathymetric highs include terrane with oceanic crustal properties, such as oceanic plateaus, seamounts, and submarine ridges, and terrane with continent crustal properties, such as microcontinents and continental fragments. Several studies on subduction with bathymetric highs have been conducted, allowing scientists to better understand its dynamic effects: 1) the transition from flat or low angle subduction to steep-slab subduction (Van Hunen et al., 2002, Van Hunen et al., 2004; Martinod et al., 2005; Gerya et al., 2009; Huangfu et al., 2016; Manea et al., 2017; Leng and Huang, 2018; Dai et al., 2020; Yan et al., 2020, 2021), 2) indentation of the trench (Dominguez et al., 1998, 2000; Morra et al., 2006), 3) surface topography development (Li F. C. et al., 2013; Ruh, 2016) and crustal shortening (Liao et al., 2018), 4) trench migration (Li Z. H. et al., 2013; Yoshida, 2017; Tao et al., 2020), 5) continental underplating (Vogt and Gerya, 2014; Magni et al., 2017). In particular, Tetreault and Buiter (2012, 2014), Yang et al. (2018) and Liu et al. (2021) have systematically presented how various crustal properties of bathymetric highs (with continent versus oceanic crustal affinity) impact the amount of accreted/subducted crust, the distinct modes of terrane accretion/complete subduction, and the deformation type of the overriding plate. However, these recent studies have employed a fixed terrane size and did not investigate the effects of terrane thickness and width on slab subduction processes.

The South China Sea is one of the largest marginal seas of the western Pacific (Deng et al., 2020). The Manila trench is located at the eastern boundary of the SCS. It was created by the subduction of the SCS plate beneath the Philippine Sea Plate (PSP) since the early Neogene, and it was induced by the northwestern plate motion of the PSP (Huang et al., 2006; Fan et al., 2016; Wu et al., 2016) (Figure 1A). For the northern Manila trench at around 18° N—21.5° N, where an obvious seaward convex is found, the morphology of the subducted SCS plate is characterized by dramatic changes from a horizontal subducting angle to near-vertical (Figure 1B) (Fan et al., 2016; Wu et al., 2016; Chen et al., 2021). A large buoyant plateau (≥300 km in width) was proposed to explain the sharp convex in the trench line, and the crustal property of the plateau was recently proposed as highly thinned continental crust (12–15 km in thickness) (Eakin et al., 2014; Lester et al., 2014; Sibuet et al., 2016; Liu et al., 2018) rather than oceanic crust as early identified (Hsu et al., 2004). For the southern Manila trench at around 15° N—18°N, the most remarkable morphological relief is an almost straight trench line and widely distributed small seamounts (≤100 km in width), such as the subducting Huangyan-Zhenbei seamounts chain. Based on high-resolution P-wave tomographic images, the morphology of the subducted SCS plate beneath the southern segment of the Manila trench showed no abrupt change in dipping angle (Figure 1C), whereas slab break-off might occur at the depths between 60 and 190 km near 17°N (Cheng et al., 2019). The trench shape and subducting slab morphology of the northern Manila trench are distinct from the southern segment. However, there have been few comprehensive analyses to discuss the relationship between the subducting bathymetric highs and the along-strike variations of the Manila trench.

(A): Tectonic setting of the Manila trench, South China Sea. Continent-ocean boundary (COB) lines indicate COB locations (modified from Liu et al., 2018). (B) and (C): Vertical cross sections of P-wave tomography along the profiles at 20°N and 17°N, respectively. The buoyant plateau and its subducted part (the gray shaded area), and the tomography sections are modified from Fan et al. (2016).

In this study, we have undertaken a series of 2-D numerical thermomechanical experiments to investigate 1) the potential crustal properties of the northeastern SCS near the Manila trench at around 20° N and 2) the key physical parameters controlling the subducting slab morphology variations between the northern and southern Manila trench. Based on our systematic numerical results, we demonstrate that dramatic changes in the dip angle in the northern Manila trench and the convex shape were caused by subduction of a large thin continental terrane, whereas the smooth morphology of subducting slab in the southern segment and straight trench were associated with normal oceanic subduction with small seamounts.

 Geological Background

The South China Sea is located in an important geodynamic intersection zone surrounded by Eurasian plates and the Pacific and Indian oceans. The partial subduction of the SCS along the Manila trench represents the last phase of a near-complete Wilson cycle, following continental rifting, breakup, and seafloor spreading. Based on deep-tow magnetic anomalies, multi-channel seismic data, the results of microfossils from IODP Expeditions and 39Ar/40Ar data, the SCS has undergone multiphase rifting events since the Late Cretaceous to Paleogene (Sun et al., 2009; Franke et al., 2014; Li et al., 2015; Sibuet et al., 2016; Ding et al., 2020; Zhang et al., 2020; Zhao et al., 2021), leading to the opening of the SCS basin at ∼32–33 Ma, and stopped spreading at ∼15 Ma in the east subbasin and ∼16 Ma in the southwest subbasin, followed by eastward subduction under the Philippine Sea Plate (PSP) along the Manila trench (Li C. F. et al., 2013, 2015; Chen et al., 2017, 2021; Jian et al., 2018; Sun et al., 2018; Deng et al., 2020; Hung et al., 2020). Large amounts of magmatism persisted for nearly 10 Ma after the cessation of seafloor spreading and generated the Zhenbei-Huangyan seamount chain (Sibuet et al., 2016; Hung et al., 2020).

There is a wide discussion on the geometry of the Manila trench (Yang et al., 1996; Bautista et al., 2001; Eakin et al., 2014; Fan et al., 2015, 2016), which forms a broad bend to the east at around 21°–18° N, trends N-S almost in a straight line from 18° to 13°N, and swerves abruptly to the ESE at its southern terminus at 13°N. Several models have been proposed to explain the sharp bend in the trench line. For example, Bautista et al. (2001) interpreted these as a collision and subsequent partial subduction of a large buoyant plateau (Figure 1). Despite its implications to the buoyancy effect, the crustal properties of this plateau in northeastern SCS remain controversial.

Geophysical studies on the deep lithospheric structures of the northeastern SCS near the northern Manila trench contribute to defining the continent-ocean boundary (COB) locations in the SCS. The crust located west of the Manila trench and north of the COB is thinned continental crust rather than thickened oceanic crust (Liu et al., 2018). Early studies defined COB1 (Figure 1) and suggested that the northeastern SCS was composed entirely of ocean crust up to 21.5°N (Hsu et al., 2004), based on E-W trending magnetic anomalies. Controversially, later evidence from refraction and multi-channel seismic (MCS) reflection data (Wang et al., 2006; Yeh et al., 2010, 2012; Lester et al., 2013, 2014; McIntosh et al., 2013, 2014; Eakin et al., 2014; Sibuet et al., 2016) defined the COB location as being more southeastward (COB2 in Figure 1). For example, Eakin et al. (2014) showed evidence for extended to hyper-extended continental crust subducting along the Manila trench, underplated to the accretionary prism at 21.5°N. Sibuet et al. (2016) has presented several features on MCS profiles suggesting that the crust of the northeastern part of the SCS is thinned continental crust intruded by post-rift volcanism. Most recently, an E-W oriented ocean bottom seismograph wide-angle refraction profile (21° N) is constructed in the northeastern SCS in 2015 (Liu et al., 2018) to further define the COB location (COB3 in Figure 1). Abundant Moho interface was shown, and a 12–15 km thick continental crust was further identified. Additionally, dramatic changes in the dipping angle of the subducted SCS plate are revealed from the northern Manila trench (20°–21.5°N). According to seismic and tomography data in this area (Fan et al., 2016), at 20°N, the SCS plate initially subducts along the Manila trench to a ∼250 km depth at a low angle of ∼25°. Then, it changes abruptly to a higher dip angle of ∼75° to a depth of ∼500 km. Further tomography studies and slab unfolding estimated the subducted portion of the SCS slab extended 400–500 km east of the present Manila trench (Wu and Suppe, 2018).

Compared to the controversial crustal properties near the northern Manila trench, the SCS plate near the southern Manila trench (at ∼15–17°N) is a typical or thickened oceanic crust, subducting with numerous seamounts. As these seamounts are formed after SCS seafloor spreading cessation, the dating results of seamount ages range from 15–6.64 Ma (Tu et al., 1992; Yan et al., 2008; Wang et al., 2009; Zhang et al., 2020) based on petrological samples. Huangyan-Zhenbei seamounts chain, located in the center of the east subbasin, is an E-W trending chain of seamounts oriented obliquely to the surrounding N055°seafloor expansion trends (Zhao et al., 2018) (Figure 1). It consists of Zhenbei seamount in the westmost (9.1 ± 0.2–10.0 ± 1.8 Ma, basalts isotopic age) (Wang et al., 2009), the Huangyan seamount in the middle, and a further east one with NEE-trending that reaches the Manila trench (Cheng et al., 2019). The crustal thickness beneath the Huangyan and Zhenbei seamounts is generally between 12 and 13.2 km based on wide-angle seismic refraction data (He et al., 2016; Zhao et al., 2018). High-resolution regional tomographic studies show that the SCS slab subducts along the southern Manila trench (at 16–17°N) at an angle of ∼45°–300 km depth, and the slab morphology is smooth without dramatic changes in the dipping angle (Figure 1C).

The 1,200 km long NS-trending Luzon island arc, stretching from 24°N to 13°N, resulted from the subduction of the SCS plate beneath the PSP along the Manila trench (Fan et al., 2016; Liu et al., 2021), which was induced by continuous NNW movement of the PSP since 25 Ma. The average plate convergence rate between the PSP and SCS was given as 7 cm/yr (Seno et al., 1993; Kreemer et al., 2003; Simons et al., 2007; Hsu et al., 2012). The Moho discontinuity is identified at a depth of 18–34 km in Luzon inferred from receiver function analysis (Besana et al., 1995). The results agree with a more recent gravity model by Manalo et al. (2015), reflecting a ∼21–31 km thick crust across the Central Philippines.

 Numerical Implementation and Model Setup Governing Equations

The momentum, continuity, and heat conservation equations for the two-dimensional creeping-flow, accounting for thermal and chemical buoyancy, were solved using modified I2VIS code (Gerya & Yuen, 2003, 2007). The incompressible continuity equation approximated conservation of mass: ∂ v x ∂ x + ∂ v z ∂ z = 0

Two-dimensional Stokes equations: ∂ σ ′ x x ∂ x + ∂ σ ′ x z ∂ z =   ∂ P ∂ x ∂ σ ′ z z ∂ z + ∂ σ ′ x z ∂ x = ∂ P ∂ z − g ρ ( T , P , C ) and a heat conservation equation: ρ C p ( D T D t ) = − ∂ q x ∂ x − ∂ q z ∂ z + H r + H a + H s q x = − k ( T , P , C ) ∂ T ∂ x ,   q z = − = − k ( T , P ,   C ) ∂ T ∂ z H r = c o n s t a n t   ( C ) , H a = T α D P D t ,   H s = σ ′ x x ε ˙ x x + σ ′ z z ε ˙ z z + 2 σ ′ x z ε ˙ x z

are used, where D/Dt is the substantive time derivative ;   k ( T , P , C ) is the thermal conductivity as a function of temperature (T), pressure (P), and composition (C) (Hofmeister, 1999); Cp is the effective isobaric heat capacity, incorporating latent heat; Hr, Ha, and HS denote radioactive heat production, the energetic effect of isothermal (de)compression (i.e., adiabatic heating/cooling), shear heating, α is the thermal expansion coefficient, σ ′ x x , σ ′ z z ,   σ ′ x z are deviatoric stress components and ε ˙ x x ,   ε ˙ z z , ε ˙ x z are strain rate components.

 Rock Rheology Implementation

Viscosity, dependence on strain rate, pressure, and temperature were defined in terms of deformation invariants: η d i f f u s i o n = A D 2 σ c r n − 1 e x p ( E a + V a P R T ) η d i s l o c a t i o n = 1 2 ( ℰ ˙ I I ) 1 − n n ( A D ) − 1 n e x p ( E a + V a P n R T ) 1 η d u c t i l e = 1 η d i f f u s i o n + 1 η d i s l o c a t i o n where ℰ ˙ I I is the second invariant of the strain rate; σ c r   denotes critical stress between dislocation creep and diffusion; A D ,   E a ,   V a ,   and   n are material constant, activation energy, activation volume, and stress exponent, respectively (Gerya, 2019; Tang et al., 2020). These material properties were determined from laboratory flow experiments and are provided in Table 1.

Material properties used in the numerical experiments.

Material State ρ0 (kg m-3) k (W m−1 K−1) Hr (μW m−3) Plastic Sin( ϕ e f f ) Viscous Flow law η0 (pan s) Ea (kJ mol−1) Va (J bar−1) n
CUC Solid 2,700 0.64 + 807/(T + 77) 1.0 0.15 Wet Quartzite 1.97 × 1017 154 1.2 2.3
CLC Solid 2,800 0.64 + 807/(T + 77) 0.25 0.15 Plagioclase_An75 4.8 × 1022 238 0.8 3.2
OUC Solid 3,000 1.18 + 474/(T + 77) 0.25 0 Wet Quartzite 1.97 × 1017 154 0.8 2.3
OLC Solid 3,000 1.18 + 474/(T + 77) 0.25 0.6 Plagioclase_An75 4.8 × 1022 238 0.8 3.2
Mantle Dry 3,300 0.73 + 1293/(T + 77) 0.022 0.6 Dry_olivine 3.98 × 1016 532 0.8 3.5
470 0.8 4.0
Wet 3,200 0 Wet_olivine 5.01 × 1020
References — 1,2 3 1 — 4 4 4 4 4

References: 1, Turcotte and Schubert, 2002; 2, Bittner and Schmeling, 1995; 3, Clauser and Huenges, 1995; 4, Ranalli, 1995. CUC, continental upper crust; CLC, continental lower crust; OUC, oceanic upper crust; OLC, oceanic lower crust.

Plasticity was implemented using the Druker-Prager yield criterion (Ranalli, 1995). The calculated creep viscosity is therefore limited as follows: σ y i e l d = C + P ⁡ sin ( φ e f f ) η p l a s t i c = σ y i e l d 2 ε ˙ I I

σ y i e l d is the yield stress. C   is the residual rock strength. P is the dynamic pressure. φ   is the internal frictional angle. φ e f f can be illustrated as the effective internal frictional angle.

With the η p l a s t i c   and η d u c t i l e , visco-plastic rheology is employed to the model where the rheology behavior depends on the minimum viscosity (Ranalli, 1995): 1 η e f f = 1 η d u c t i l e + 1 η p l a s t i c

Model Setup

The initial configuration of reference 2D numerical model is shown in Figure 2. The numerical model box with 1,361 by 351 nodal points is non-uniform and corresponded to 4,000 by 1,400 km physical dimension. The rectangular grid contain a 1,000 km wide high-resolution area of 1 by 1 km grid step size in the center of the domain. The rest of the model is at a lower resolution (up to 10 by10 km grid step size). Over 5.7 million Lagrangian markers randomly distribute in the whole model domain.

Initial setups of the numerical models (A) Full box initial geometry (4,000 × 1,400 km). The orange arrow denotes the right overthrusting velocity. (B) Zoomed area of the collision domain (900 × 200 km). The white solid lines are isotherms with an increment from 100 to 1,300°C. Composition color codes: 0, stick air; 1, water; 2, sediments; 3 and 4, upper and lower continental crust, respectively; 5 and 6, upper and lower oceanic crust, respectively; 7, lithospheric mantle; 8, asthenosphere mantle; 9. hydrated mantle (weak zone). (C) and (D) Lithosphere yield stress profiles for continental terrane and oceanic plateau, respectively. See Table 1 for the rheological parameters. UC, upper crust; LC, lower crust.

The oceanic crust contain a thin continental terrane, fated to collide with the seaward-moving overriding continental plate. The oceanic crust is composed of 2-km-thick upper crust of hydrothermally-altered basalt, underlain by 5-km-thick lower crust of gabbroic rocks that covered 2,500 km horizontally. The continental crust is felsic and has a total thickness of 30 km, composed of 15 km upper and 15 km lower crust that extend over 1,500 km. The large 300 km wide and 15 km thick continental terrane is defined 200 km from the trench on the oceanic plate. Oceanic plateau (Figure 2B) paired numerical experiments with similar parameters differing only by the crustal property of the terrane are run in parallel to compare the different dynamic effects with continental terrane. Lithosphere yield stress profiles comparing the properties of continental terrane and oceanic plateau are shown in Figures 2C,D. Both the asthenosphere and upper mantle are composed of anhydrous peridotite and are defined by the temperature profile. The rheological parameters used in the reference model are summarized in Table 1. We keep our model with a far-field push on the overriding plate (7 cm/yr) to be similar to the average plate convergence rates of the Manila trench (Seno et al., 1993; Kreemer et al., 2003; Simons et al., 2007; Hsu et al., 2012) while simultaneously as simple as possible to test different model parameters. It is also important to note that prescribing velocity of overlying rather than subducting plate enhances tendency of slab flattening (van Hunen et al., 2004).

All mechanical boundary conditions are free slip, and only the lower boundary was permeable, satisfying an external free slip boundary condition (Gorczyk et al., 2007; Ueda et al., 2008). In addition, on the top of the rocky portion of the model is treated as an internal free surface (Schmeling et al., 2008) by using a top layer (of 20 km thickness) with low viscosity (1018 Pas) and low density (1 kg/m3 for air, 1,000 kg/m3 for sea water) to allow for the surface topographic evolution.

The initial temperature field of the oceanic plate is defined by its oceanic geotherm for a specific lithospheric cooling age of 20 Ma, according to the subducting slab age of the Manila trench (Wu et al., 2016). The oceanic plateau embedded into the oceanic crust is assumed to have the same thermal structure as the oceanic lithosphere. Therefore, the initial temperature field of the continental plate is increased linearly from 0°C at the surface to 1,344 °C at the lithosphere-asthenosphere boundary (140 km depth). For the asthenospheric mantle (>140 km), a thermal gradient of 0.5°C km−1 is used.

 Subduction of Thin Continental Terrane

A series of models (Model-C100, Model-C200, Model-C300, Model-C400) are conducted to test the influence of a bathymetric high with continental affinity. We focus on the thickness and width of the terrane. Four groups of experiments with fixed terrane width (100 km, 200 km, 300 km, and 400 km for each group) are examined by changing the terrane thickness (10 km, 15 km, 20 km, 25 km, and 30 km in thickness for each group). A representative selection of the models is shown in Table 2 to discuss how different size terranes embedded in the subducting slab affect the subduction process.

Parameters and results of representative experiments.

Model name Terrane Terrane Width (km) Terrane Thickness (km) Results
Model-C300-10 Continent (large) 300 10 Figure 4
Model-C300-15 (Reference model) Continent (large) 300 15 Figure 3
Model-C300-20 Continent (large) 300 20 Figure 4
Model-C300-25 Continent (large) 300 25 Figure 4
Model-C300-30 Continent (large) 300 30 Figure 4
Model-C100-10 Continent (small) 100 10 Figure 6
Model-C100-15 Continent (small) 100 15 Figure 6
Model-C100-20 Continent (small) 100 20 Figure 6
Model-C100-25 Continent (small) 100 25 Figure 6
Model-C100-30 Continent (small) 100 30 Figure 6
Model-O300-15 Oceanic (large) 300 15 Figure 7
Model-O100-10 Oceanic (small) 100 10 Figure 8
Model-O100-15 Oceanic (small) 100 15 Figure 8
Model-O100-20 Oceanic (small) 100 20 Figure 8
Model-O100-25 Oceanic (small) 100 25 Figure 8
Model-O100-30 Oceanic (small) 100 30 Figure 8
Model-O400-10 Oceanic (large) 300 10 Figure 8
Model-O400-15 Oceanic (large) 300 15 Figure 8
Model-O400-20 Oceanic (large) 300 20 Figure 8
Model-O400-25 Oceanic (large) 300 25 Figure 8
Model-O400-30 Oceanic (large) 300 30 Figure 8
Reference Model

The numerical evolution of the reference model (Model-C300-15) is shown in Figure 3. In this model, a thin continental terrane (15-km thick, 300-km wide) is embedded in the oceanic plate. The slab initially subducts with a smooth morphology, and the thin continental terrane passes through the subduction channel (Figure 3A). A small part of the upper crust of the terrane is accreted to the overriding continent margin. This leads to an uplift in the accretionary wedge (Figure 3B1). Then the terrane is subducted into the deep mantle creating a shallower slab angle at a depth of ∼100 km. The downgoing oceanic plate is substantially narrowed at a depth of ∼250 km and shows a rheologically weak gap in the narrowed part (Figure 3C2), leading to an abrupt change in the dipping angle (from ∼20° to ∼58°) (Figure 3C1). Finally, the slab pull causes the occurrence of slab break-off (Figure 3D). This model illustrates that the subduction of a large (300-km wide) terrane with thin continental (15-km thick) affinity leads to an abrupt change in the dip angle and an uplift in the accretionary prism.

Evolution of the reference model (Model-C300-15), with thinned continental terrane (15-km thick, 300-km wide). Left: Evolution of composition in Model-C300-15. The composition code shown here is the same as in Figure 1. Right: Evolution of viscosity in Model-C300-15. Black arrows in (A1), (B1), (D1) and (C2) indicates the partially subducted continental terrane, the uplift in the accretionary wedge, the slab break off, and the rheologically weak gap in the mantle, respectively. Red box in (B1) indicates the uplift in the accretionary wedge. The magnitude of angle in Figure C1 shows the abrupt change of dipping angle.

 Variations in Thickness of Large-Size Terranes (≥300 km in Width)

Model-C400 differs from Model-C300 only because the continental terrane has a larger width (400-km wide continental terrane). The evolution of Model-C400 is similar to Model-C300. Figure 4 shows the morphology of the subducted plate right before breaking off.

Comparison of large-size terrane models with varied terrane thickness, showing the slab morphology right before breaking off. (A–E) Slab morphology of Model-C300, in which the terrane is 300 km wide. The terrane thickness is 10–30 km from A to E, respectively. (C1) Topography result of Figure (C) (F–J) Slab morphology of Model-C400, in which the terrane is 400 km wide. The terrane thickness is 10–30 km from F to J, respectively. (H1) Topography result of Figure (H) and The black arrows in (A), (B), (F), and (G) indicate the change in dipping angle. The black arrows in C and H refer to no abrupt change in dipping anlge before slab break off. The red boxes and arrows indicate the uplift in the accretionary wedge. The black arrows in (D), (E), (I), and (J) show the blocked terranes.

For terranes with thin continent crust (≤15 km in thickness) (Liu et al., 2018; Deng et al., 2020), abrupt changes in dipping angle are shown in Figures 4A,B,F,G, and the plates narrow at a depth of 200–300 km. Figures 4A,B show that the significant change in the dip angle happens after the terrane pushes into the mantle in Model-C300, whereas part of the terrane is still going through the subduction channel in Model-C400 due to a larger width.

In the case of a normal-thickness continental terrane (25–30 km in thickness), most of the terrane is blocked in the subduction zone before the slab breaks off, and only smaller parts are sheared off and dragged down into the mantle (Figures 4D,E,I,J). For both Model-C300 and Model-C400, the oceanic slab subducts at a shallow angle at an early stage and gradually steepens to nearly vertical. The morphology of the downgoing plate remained smooth and showed no abrupt change in the dipping angle before the slab broke off. Most parts of the thick continental terrane were blocked due to its buoyancy (Cloos, 1993), and terrane collision and lateral accretion occurred to accommodate the constant convergence (Figure 5). A “subduction zone jump” (Yan et al., 2021) after the collision is shown in Figure 5, which may be caused by the detachment of the buoyant crust of the terrane (Zhang et al., 2020).

Evolution of Model-C300-30, with normal-thickness continental terrane (30-km thick, 300-km wide). Left: Evolution of composition in Model-C300-30. The composition code shown here is the same as in Figure 1. Right: Evolution of viscosity in Model-C300-30. Black arrows in (A1), (B1), and (C1) indicate the slab break-off, the blocked terrane and the initiation of a new subduction zone behind it, and the location of the new subduction zone, respectively.

Terranes with thin continent crust (20 km in thickness) act as a transition from thin terrane subduction with an abrupt morphology change to normal-thickness terrane collision. The oceanic plate is subducted at a gradually deeper angle, and there is no abrupt change in slab morphology before the slab breaks off (Figures 4C,H). Most of the upper crust of the subducting thin continental terrane is accreted to the leading edge of the overriding plate. In contrast, the lower crust and minor part of the upper one bypass the accretionary prism and is lost by subduction. As a result, crustal thickening accumulate in the downgoing plate and propagate toward the accretionary wedge. The depth at which part of the terrane detached from the subducting slab is as shallow as 30 km. Thickened upper crust contributes to a vertical increase of accretionary prism (Figures 4C1,H1), which is similar to the uplift in the reference model (Figure 3B).

The Model-C300 and Model-C400 sets demonstrate that the thickness of a large continental terrane plays a vital role in slab morphology during the subduction process. A thin continental terrane causes an abrupt change in dipping angle before slab breakoff, while a normal-thickness continental terrane resists subduction and encourages “subduction zone jump” (Figure 5).

 Variations in Thickness of Moderate/Small-Size Terranes (200 km/100 km in Width)

Model-C100 and Model-C200 are identical to the reference model except for the terrane width. Model-C100 has a small-size continental terrane (100-km wide) (Yan et al., 2021), and Model-C200 has a moderate-size continental terrane (200-km wide) (Gerya et al., 2009). Figure 6 shows the morphology of the subducting plate right before breaking off.

Compared small/moderate-size continental terrane models with varied terrane thickness, showing the slab morphology right before breaking off. (A–E) Slab morphology of Model-C100, in which the terrane is 100 km wide. The terrane thickness is 10–30 km from A to E, respectively. (F–J) Slab morphology of Model-C200, in which the terrane is 200 km wide. The terrane thickness is 10–30 km from F to J, respectively. The black arrows in (A), (B), (F) and (G) indicate the change in dipping angle and the subducted terrane. The black arrows in (C), (D), and (H) indicate the partially accreted terrane crust. The red boxes in (C), (D), and (H) indicate the uplift in accretionary wedge. The black arrows in (E), (I), and (J) indicate the blocked terranes.

A small and thin continental terrane (100 km in width and ≤ 15 km in thickness) leads to complete terrane subduction. As the whole terrane is dragged down into the deep mantle, the downgoing plate show no change in the dipping angle (Figure 6A) or a relatively slight change from 29° to 44° (Figure 6B). The change in slab morphology is more obvious in the moderate-size continental terrane subducting process (Figures 6F,G). When the small/moderate-size terrane has a thicker crust (20–30 km thick), the model evolution is similar to Model-C300 and Model-C400, which present the large continental terrane subduction. The subducting slab show no abrupt change in the dipping angle before breaking off (Figures 6C–F, H–J), with blocked terrane and subsequent “subduction zone jump” in normal-thickness continental terrane models (25–30 km thick), and partial terrane accretion in thin continental terrane model (20 km thick).

Model-C100 and Model-C200 illustrate that the width of the continental terrane is also a controlling factor on the subducting slab morphology. Thus, the moderate-size continental terrane or wider ones are more likely to form the abrupt change of dipping angle during the subducting process.

 Subduction of Oceanic Plateau

We test the effect of oceanic plateau subduction to clarify its different influence on slab morphology from continental terrane subduction. Figure 7 show the evolution of Model-O300-15 with an oceanic plateau (15-km thick, 300-km wide) embedded in the subducting plate. The oceanic plate begin to subduct along the weak zone. Once the oceanic plateau enter the subduction zone, a large part of the upper crust of the plateau is sheared off and accreted into the accretionary prism (Figures 7A1,B1). When the oceanic plateau bypass the subduction channel, the downgoing oceanic plate show no abrupt change in dipping angle (Figures 7A1, B1, C1). Materials of the forearc region of the overriding plate do not show obvious uplift. The slab break-off occur before the whole plateau is dragged down into the subduction channel (Figure 7C1). Significant flexural stresses are generated in the subducting plate before breaking off (Figure 7B2).

Evolution of Model-O300-15, with oceanic plateau (15-km thick, 300-km wide). Left: Evolution of composition in Model-O300-15. The composition code shown here is the same as in Figure 1. Right: Evolution of viscosity in Model-O300-15. Black arrows in (A1), (B1), and (C1) indicate the subducting oceanic plateau, the partially accreted plateau crust and the smooth subducting slab morphology, and the slab break-off, respectively.

A series of oceanic plateau models (Model-O100, Model-O200, Model-O300, and Model-O400) are performed to compare their influence on subducting slab morphology and corresponding thicknesses to continental terrane models (Model-C100, Model-C200, Model-C300, and Model-C400). Here a selection of the models performed is shown in Figure 8 to compare models with small/large-size oceanic plateaus (100 km/400 km in width). In Figure 8, it is clear that neither a small oceanic plateau nor a large one favored the abrupt change in the dipping angle, irrespective of how thick the oceanic plateau is. For plateau thicknesses no more than 20 km, part of the upper crust is accreted into the accretionary wedge, and the rest of the plateaus subduct into the deep mantle (Figures 8A–C,F–H). For plateau with an over-thickened crust (25 km or 30 km in thickness), most part of the plateau is blocked in the subduction zone before slab break-off (Figures 8D,E,I,J). According to our model results, a thicker oceanic plateau on the subducting plate favor a faster break-off event.

Compared small/large-size oceanic plateau models with varied plateau thickness, showing the slab morphology right before breaking off. (A–E) Slab morphology of Model-O100, in which the oceanic plateau is 100 km wide. The plateau thickness is 10–30 km from A to E, respectively. (F–J) Slab morphology of Model-O400, in which the oceanic plateau is 400 km wide. The plateau thickness is 10–30 km from F to J, respectively. All the model results in this figure show no abrupt change in dipping angle.

These model results indicate that the oceanic plateau subduction is more likely to form a smooth slab morphology, differing from the abrupt change in the dipping angle in thin continental terrane subduction models.

 Discussion The Role of Crustal Properties and Size of Bathymetric Highs in Subducting Slab Morphology

Figure 9 summarizes the model results from varied terrane crustal properties, width, and thickness. The reference model has a 300 km wide and 15 km thick continental terrane. It is characterized by deep continental crust subduction and by an abrupt change in subducting slab morphology (Δθ > 20°, Figure 9). Additional tests show that: 1) terrane width or thickness variation has no significant effect on subducting slab morphology of oceanic plateau models (Model-O100, Model-O200, Model-O300, and Model-O400). No change in the dipping angle is observed (Δθ = 0°) before slab break-off occurs (Figures 8, 9). 2) continental terrane models with a terrane thickness ≥20 km (Figures 4C–E, H–J, Figures 6C–E, H–J and Figure 9) showed no abrupt change in the dipping angle before shallow slab break off, whatever how wide the terrane is. 3) in the context of thin continental terrane models (terrane thickness = 10 km or 15 km) (Figures 4A,B, F–G), Figure 6A,B, F–G and Figure 9), large-size terrane (Model-C300 and Model-C400) favors for the abrupt change of dipping angle (Δθ > 20°, Figure 9) during the subducting process. In contrast, small-size terrane (Model-C100) is more likely to result in a slight change in dipping angle (Δθ < 20°, Figure 9) or no change (Δθ = 0°, Figure 9) before slab break-off. Moderate-size terrane (Model-C200) is a transition type from abrupt change to slight change, and its Δθ = 20° (Figure 9).

Regime diagram showing the effect of terrane thickness and width on slab morphology. Each colored solid line and dot indicate one set of model with a fixed width (e.g., the green line, model C-100, represents continental terrane model with a terrane width of 100 km; the purple hollow circle, model O-200, represents oceanic plateau model with a terrane width of 200 km). It is worth noting that oceanic plateau models (O-100, O-200, O-300, O-400) are all with Δθ = 0°, so only colored dots are shown on horizontal axis rather than lines. Three distinct modes of subducting slab morphology are observed: abrupt change in dipping angle (Δθ > 20°), slight change in dipping angle (Δθ < 20°), no change in dipping angle before slab break-off (Δθ = 0°). The horizontal axis depends on terrane thickness. The vertical axis depends on Δθ. Δθ = θ2—θ1. “C” and “O” represent continental terrane model and oceanic plateau model, respectively.

Our oceanic plateau models show no abrupt change in the dipping angle before slab break-off, consistent with the results of Cheng et al. (2019), who showed that seamount subduction promotes the break-off process, and Gerya et al. (2009), who showed smooth morphology of slab position lines. Other previous numerical models with bathymetric highs, including continental fragments, seamounts, oceanic plateaus, and island arcs, focus on various types of accretion (Tetreault and Buiter, 2012, 2014; Li Z. H. et al., 2013; Vogt and Gerya, 2014; Yang et al., 2018). This partially compares well with our numerical experiments. Part of the upper crust of the subducting thin continental terrane/oceanic plateau is accreted to the leading edge of the overriding plate. The above models are 20 km thick and 100–200 km wide (small to moderate-size) bathymetric highs, consistent with the average thickness of global thin continental fragments/oceanic plateaus. Recent numerical modeling studies test small to large-size oceanic plateaus (Tao et al., 2020; Yan et al., 2021; Almeida et al., 2022; Wang et al., 2022) to investigate how buoyant plateaus contribute to dip angle change, subduction polarity reversal, and subduction zone jumping. Previous analogue models investigate the effects of seamount subduction on the structural deformation of the accretionary wedge, especially the evolution of faults in the wedge (Dominguez et al., 1998, 2000; Li F. C. et al., 2013; Wang et al., 2021). Li F. C. et al. (2013) investigates the seamount subduction along the Manila trench by combining analog and gravity modeling, and suggests that the dip anlge of subducting plate may be associated with extict mid-ocean ridge subduction. However, as most of the analogue models focus on seamount or oceanic ridge subduction, the effects of continental terrane subduction are rarely to be compared with in analoge modeling. Our models systematically investigate the subduction process of different bathymetric highs with varied crustal properties, width, and thickness and further illustrate that the abrupt change of subducting slab morphology is not only controlled by the crustal properties of bathymetric high but also by its thickness and width.

Sharp changes in subducting slab angle imply localized slab bending, which is driven by combined effects of rheological weakness and positive buoyancy of subducted continental crust. Similar localized bending process (segmentation) of purely oceanic subducting slabs has been recently suggested by Gerya et al. (2021), which is primarily driven by strain weakening of outer-rise normal faults, and grain-size reduction (ductile damage) of the lithospheric mantle. However, our models do not consider the grain-size reduction in the slab, which prevents us from directly comparing our models with seismic tomography data beneath the Manila trench. Such comparison also requires (cf. discussion in Gerya et al., 2021) better resolution of seismic data in Manila subduction zone, which is currently unavailable.

 Implications for the Subduction Along the Manila Trench

The Manila trench is characterized by distinct trench shape and subducting slab morphology between the northern and southern segment. Here, we apply our results to understand the relationship between the subducting bathymetric highs and the along-strike variations of the Manila trench. This study performed an overriding plate push at the rate of 7 cm/yr (Seno et al., 1993; Kreemer et al., 2003; Simons et al., 2007; Hsu et al., 2012), a relatively young subducting oceanic plate (20 Myrs) (Wu et al., 2016), and a thin continental terrane (15 km thick in the reference model) (Wang et al., 2006; Yeh et al., 2010, 2012; Lester et al., 2013, 2014; McIntosh et al., 2013, 2014; Eakin et al., 2014; Sibuet et al., 2016; Liu et al., 2018) to compare with the geological settings of the Manila trench.

As shown in Figure 1, at around 20°N, where the maximum seaward convex of the Manila trench is found, the crustal property of the northeastern SCS was initially characterized as the oceanic crust (Hsu et al., 2004), whereas later studies suggested that it is thin continent crust with a thickness of ∼12–15 km (Wang et al., 2006; Yeh et al., 2010, 2012; Lester et al., 2013, 2014; McIntosh et al., 2013, 2014; Eakin et al., 2014; Sibuet et al., 2016; Liu et al., 2018). Additionally, dramatic changes in the dipping angle of the subducted SCS plate are revealed from the northern Manila trench (20° N) (Fan et al., 2016; Wu et al., 2016; Chen et al., 2021). According to seismic and tomography data in this area (Fan et al., 2016), at 20° N, the SCS plate subducts initially along the Manila trench to ∼250 km depth at a low angle of ∼25°. Then, it abruptly changes to a higher dip angle of ∼75° to a depth of ∼500 km (Figure 10A2). In our model results, the thin continental terrane was initially subducted along the subducting channel at a low angle. Then, it changed abruptly to a much higher dip angle (Figure 10A1), whereas all oceanic plateau models do not result in an abrupt change in the dipping angle. Thus, our continental terrane model result is consistent with the tomographic profile at the northern Manila trench (20°N), and further support that the crust located west of the Manila trench and around 20°N is a thin continental crust, rather than oceanic plateau.

Comparison of modeled and observed subducting South China Sea (SCS) slab morphology along the Manila trench. (A1) Model results of thin continental terrane subduction. The red arrow indicates an abrupt change in the dip angle. (A2) Vertical cross section of P-wave tomography along the profile at northern Manila trench (20° N). Modified from Fan et al., 2016. (A3) Topography evolution of thin continental terrane subduction. (B1) Model results of oceanic plateau subduction. (B2) Vertical cross section of P-wave tomography along the profile at southern Manila trench (17° N). Modified from Fan et al., 2016. (B3) Topography evolution of oceanic plateau subduction. The red circle indicates the vertical increase in the accretionary prism volume.

In contrast, the southern Manila trench is characterized by widely distributed small-size seamounts (<100 km in width), e.g., Huangyan-Zhenbei seamounts chain with a crustal thickness between 13 and 14 km (He et al., 2016; Zhao et al., 2018; Cheng et al., 2019). Thus, our model with a 15 km-thick and 100 km-wide oceanic plateau is suitable for this area. This model does not result in an abrupt change in dipping angle (Figure 10B1), which is also compared well with the smooth slab morphology in the tomographic profile of the southern Manila trench (17° N) (Figure 10B2). Futher investigation of lateral slab morphology variations will require application of 3D thermomechanical modelling approaches. Also, the model results indicate that the most important factor controlling the occurrence of an abrupt change in dipping angle is the crustal properties of bathymetric high along the Manila trench. Only subducting plate with continental terrane may initially subduct along the trench at a low angle and then changes abruptly to a higher dip angle, while comparable size oceanic plateau subduction favors for smooth subducting slab morphology.

Our model results exhibit that the second-order factor controlling the subducting slab morphology is the thickness of continental terrane. A thin continental terrane (10 km or 15 km thick) has a strong ability to result in an abrupt change in dipping angle before slab breakoff (Figures 4, 6, 9). In comparison, a thicker continental terrane (≥25 km thick) is more likely to be blocked in the subduction channel (Figures 4, 6), and favors for “subduction zone jump.” The latter point is similar to the ‘subduction zone jump’ phenomenon in Tao et al. (2020) and Yan et al. (2021), which is beyond the focus of this paper. In addition, the continental terrane width also has implications for affecting the downgoing slab morphology. According to our model results, the moderate-size continental terrane or wider ones (≥200 km in width) are more likely to form the abrupt change of dipping angle during the subducting process (Figures 4, 6, 9). These results are comparable with geological observations of the Manila trench, where large-size (≥300 km in width) and thin (≤15 km in thickness) continental terrane subduction leads to an abrupt change in the dipping angle in the northern segment, and small-size (≤100 km in width) seamounts subduction results in a smooth subducting slab morphology in the southern segment.

The geophysical studies of northern Manila trench indicate that the accretionary wedge shows an vertical volume increase, and it may be due to the accreted crustal material from the thin continental terrane (Eakin et al., 2014). According to our models, Figures 10A3, B3 compare the topography evolution of large thin continental terrane subduction with small oceanic plateau subduction. The former contributes to a relatively evident vertical increase in the accretionary wedge volume (Figure 10A3), which is similar with the topographic uplift in Figures 4C1, H1. These topographic results of continental terrane models agree with the geophysical observations from Eakin et al. (2014).

 Conclusion

We test a series of 2-D geodynamic models to investigate subduction processes involving thin continental terrane and oceanic plateau and provide insights from numerical modeling on properties of the subducted crust of the South China Sea along the Manila trench. There are three key findings of the study.

1. The first-order factor controlling an abrupt change in the dipping angle is the crustal properties of bathymetric highs. Subducting plate with continental terrane initially subducts along the trench at a low angle and then changes abruptly to a higher dip angle. In contrast, comparable size oceanic plateau subduction does not result in the abrupt change in the dipping angle.

2. For continental terranes, the crustal thickness and terrane width affect the subducting slab morphology. The subduction of a wide continental terrane (≥300 km) with thin crust (≤15 km in thickness) favors the abrupt change in dipping angle. Overthickened terranes (≥25 km in thickness) are more likely to be blocked in the subduction zone.

3. The model results explain the differences in subducting slab morphology between the northern (around 20° N) and southern (around 17°N) segments of the Manila trench. For the northern Manila trench, numerical models with large thin continental terrane (≥300 km in width and ≤15 km in thickness) lead to an abrupt change in dipping angel, which corresponds to tomography profile at 20°N. On the other hand, for the southern Manila trench, models with small (≤100 km in width) oceanic plateaus (seamounts) are characterized by smooth subducting slab morphology, which corresponds to the tomography profile at 17°N.

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

LM and JL conceived the study. LM and LC designed the numerical experiments. LM performed numerical experiments. LM, LC, ZC, and TG analyzed and interpreted results. TG provided the I2VIS code and guidance on improving the paper. LM, LC, ZC, TG, and JL. wrote the paper.

 Funding

This work is supported by Natural Science Foundation of Zhejiang, China (LQ20D060002), the Scientific Research Fund of the Second Institute of Oceanography, MNR, China (JG 1907), the Innovation Group Project of Southern Marine Science and Engineering Guangdong Laboratory (Zhuhai) (No. 311020018).

 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.

The authors would like to thank Dr. Jie Liao for improving the model setup and commenting the early draft. The authors would like to thank Dr. Fucheng Li and one reviewer for their insightful and detailed reviews which substantially improved the manuscript. The authors also would like to thank Shangguo Chen for improving Figure 2.

 References Almeida J. Riel N. Rosas F. M. Duarte J. C. Kaus B. (2022). Self-replicating Subduction Zone Initiation by Polarity Reversal. Commun. Earth Environ. 3, 55. 10.1038/s43247-022-00380-2 Bautista B. C. Bautista M. L. P. Oike K. Wu F. T. Punongbayan R. S. (2001). A New Insight on the Geometry of Subducting Slabs in Northern Luzon, Philippines. Tectonophysics 339, 279–310. 10.1016/S0040-1951(01)00120-2 Besana G. M. Shibutani T. Hirano N. Ando M. Bautista B. Narag I. (1995). The Shear Wave Velocity Structure of the Crust and Uppermost Mantle beneath Tagaytay, Philippines Inferred from Receiver Function Analysis. Geophys. Res. Lett. 22, 3143–3146. 10.1029/95GL03319 Bittner D. Schmeling H. (1995). Numerical Modeling of Melting Processes and Induced Diapirism in the Lower Crust. Geophys. J. Int. 123, 59–70. 10.1111/j.1365-246X.1995.tb06661.x Chen L. Hu J. W. Yang D. H. Song H. B. Wang Z. H. (2017). Kinematic Models for the Opening of the South China Sea: An Upwelling Divergent Flow Origin. J. Geodyn. 107, 20–33. 10.1016/j.jog.2017.03.002 Chen S. S. Gao R. Wang Z. W. Hou T. Liao J. Yan C. (2021). Coexistence of Hainan Plume and Stagnant Slab in the Mantle Transition Zone beneath the South China Sea Spreading Ridge: Constraints from Volcanic Glasses and Seismic Tomography. Lithosphere (Special 2), 1–22. 10.2113/2021/6619463 Cheng Z. H. Ding W. W. Faccenda M. Li J. B. Lin X. B. Ma L. T. (2019). Geodynamic Effects of Subducted Seamount at the Manila Trench: Insights from Numerical Modeling. Tectonophysics 764, 46–61. 10.1016/j.tecto.2019.05.011 Clauser C. Huenges E. (1995). “Thermal Conductivity of Rocks and Minerals,” in Rock Physics and Phase Relations: A Handbook of Physical Constants. Editor Ahrens T. J. (American Geophysical Union), 105–126. Cloos M. (1993). Lithospheric Buoyancy and Collisional Orogenesis: Subduction of Oceanic Plateaus, Continental Margins, Island Arcs, Spreading Ridges, and Seamounts. Geol. Soc. Am. Bull. 105, 715–737. 10.1130/0016-7606(1993)105<0715:lbacos>2.3.co;2 Dai L. M. Wang L. L. Lou D. Li Z. H. Dong H. Ma F. F. (2020). Slab Rollback versus Delamination: Contrasting Fates of Flat-Slab Subduction and Implications for South China Evolution in the Mesozoic. J. Geophys. Res. Solid Earth 125, e2019JB019164. 10.1029/2019JB019164 Deng H. D. Ren J. Y. Xiong P. Rey P. F. McClay K. R. Watkinson I. M. (2020). South China Sea Documents the Transition from Wide Continental Rift to Continental Break up. Nat. Communs. 11, 4583. 10.1038/s41467-020-18448-y Ding W. W. Sun Z. Mohn G. Nirrengarten M. Tugend J. Manatschal G. (2020). Laterl Evolution of the Rift-To-Drift Transition in the South China Sea: Evidence from Multi-Channel Seismic Data and IODP Expeditions 367&368 Drilling Results. Earth Planet. Sci. Lett. 351, 115932. 10.1016/j.epsl.2019.115932 Dominguez S. Lallemand S. E. Malavieille J. von Huene R. (1998). Upper Plate Deformations Associated with Seamount Subduction. Tectonophysics 293, 207–224. 10.1016/S0040-1951(98)00086-9 Dominguez S. Malavieille J. Lallemand S. E. (2000). Deformation of Accretionary Wedges in Response to Seamount Subduction: Insights from Sandbox Experiments. Tectonics 19 (1), 182–196. 10.1029/1999tc900055 Eakin D. H. McIntosh K. D. Van Avendonk H. J. A. Lavier L. Lester R. Liu C. S. (2014). Crustal-scale Seismic Profiles across the Manila Subduction Zone: The Transition from Intraoceanic Subduction to Incipient Collision. J. Geophys. Res. Solid Earth 119, 1–17. 10.1002/2013JB010395 Fan J. K. Zhao D. P. Dong D. D. (2016). Subduction of a Buoyant Plateau at the Manila Trench: Tomographic Evidence and Geodynamic Implications. Geochem. Geophys. Geosyst. 17, 571–586. 10.1002/2015GC006201 Fan J. Wu S. G. Spence G. (2015). Tomographic Evidence for a Slab Tear Induced by Fossil Ridge Subduction at Manila Trench, South China Sea. Int. Geol. Rev. 57, 998–1013. 10.1080/00206814.2014.929054 Franke D. Savva D. Pubellier M. Steuer S. Mouly B. Auxietre J. L. (2014). The Final Rifting Evolution in the South China Sea. Mar. Pet. Geol. 58, 704–720. 10.1016/j.marpetgeo.2013.11.020 Gerya T. V. Fossati D. Cantieni C. Seward D. (2009). Dynamic Effects of Aseismic Ridge Subduction: Numerical Modelling. Eur. J. Mineral. 21, 649–661. 10.1127/0935-1221/2009/0021-1931 Gerya T. V. (2019). Introduction to Numerical Geodynamic Modelling. Second Edition. London: Cambridge Univ. Press, 73–80. Gerya T. V. Yuen D. A. (2003). Rayleigh – Taylor Instabilities from Hydration and Melting Propel ‘Cold Plumes’ at Subduction Zone. Earth Planet. Sci. Lett. 212, 47–62. 10.1016/s0012-821x(03)00265-6 Gerya T. V. Yuen D. A. (2007). Robust Characteristics Method for Modeling Multiphase Visco- Elasto- Plastic Thermo-Mechanical Problems. Phys. Earth Planet Inter. 163, 83–105. 10.1016/j.pepi.2007.04.015 Gorczyk W. Willner A. P. Gerya T. V. Connolly J. A. D. Burg J. P. (2007). Physical Controls of Magmatic Productivity at Pacific-type Convergent Margins: Numerical Modelling. Phys. Earth Planet Inter. 163, 209–232. 10.1016/j.pepi.2007.05.010 He E. Y. Zhao M. H. Qiu X. L. Sibuet J. -C. Wang J. Zhang J. Z. (2016). Crustal Structure across the Post-spreading Magmatic Ridge of the East Sub-basin in the South China Sea: Tectonic Significance. J. Asian Earth Sci. 121, 139–152. 10.1016/j.jseaes.2016.03.003 Hofmeister A. M. (1999). Mantle Values of Thermal Conductivity and the Geotherm from Phonon Lifetimes. Science 283 (5408), 1699–1706. 10.1126/science.283.5408.1699 Hsu S. K. Yeh Y. C. Doo W. B. Tsai C. H. (2004). New Bathymetry and Magnetic Lineations Identifications in the Northernmost South China Sea and Their Tectonic Implications. Mar. Geophys. Res. 25, 29–44. 10.1007/s11001-005-0731-7 Hsu Y. J. Yu S. B. Song T. R. A. Bacolcol T. (2012). Plate Coupling along the Manila Subduction Zone between Taiwan and Northern Luzon. J. Asian Ear. Sci. 51, 98–108. 10.1016/j.jseaes.2012.01.005 Huang C. Y. Yuan P. B. Tsao S. J. (2006). Temporal and Spatial Records of Active Arc-Continent Collision in Taiwan: a Synthesis. Bull. Geol. Soc. Am. 118 (3-4), 274–288. 10.1130/B25527.1 Huangfu P. P. Wang Y. J. Cawood P. A. Li Z. H. Fan W. M. Gerya T. V. (2016). Thermo-mechanical Controls of Flat Subduction: Insights from Numerical Modeling. Gondwana Res. 40, 170–183. 10.1016/j.gr.2016.08.012 Hung T. D. Yang T. Le B. M. Yu Y. Q. Xue M. Liu B. H. (2020). Curstal Structure across the Extinct Mid-ocean Ridge in South China Sea from OBS Receiver Functions: Insights into the Spreading Rate and Magma Supply Prior to the Ridge Cessation. Geophys. Res. Lettetrs 48, e2020GL089755. 10.1029/2020GL089755 Jian Z. Larsen H. C. Alvarez Zarikian C. A. (2018). Expedition 368 Preliminary Report: South China Sea Rifted Margin. Texas, United States: International Ocean Discovery Program. Kreemer C. Holt W. E. Haines A. J. (2003). An Integrated Global Model of Present-Day Plate Motions and Plate Boundary Deformation. Geophys. J. R. Astron. Soc. 154 (1), 8–34. 10.1046/j.1365-246X.2003.01917.x Leng W. Huang L. Z. (2018). Progress in Numerical Modeling of Subducting Plate Dynamics. Sci. China Earth Sci. 61, 1761–1774. 10.1007/s11430-017-9275-4 Lester R. Mcintosh K. Van Avendonk H. J. A. Lavier L. Liu C. S. Wang T. K. (2013). Crustal Accretion in the Manila Trench Accretionary Wedge at the Transition from Subduction to Mountain-Building in Taiwan. Earth Planet. Sci. Lett. 375, 430–440. 10.1016/j.epsl.2013.06.007 Lester R. Van Avendonk H. J. A. Mcintosh K. Lavier L. Liu C. S. Wang T. K. (2014). Rifting and Magmatism in the Northeastern South China Sea from Wide-Angle Tomography and Seismic Reflection Imaging. J. Geophys. Res. Solid Earth 119, 2305–2323. 10.1002/2013JB010639 Li C. F. Li J. B. Ding W. W. Franke D. Yao Y. J. Shi H. S. (2015). Seismic Stratigraphy of the Central South China Sea Basin and Implications for Neotectonics. J. Geophys. Res. Solid Earth 120, 1377–1399. 10.1002/2014JB011686 Li C. F. Lin J. Kulhanek D. K. (2013c). Opening of the South China Sea and its Implications for Southeast Asian Tectonics, Climates, and Deep Mantle Processes since the Late Mesozoic. Texas, United States: International Ocean Discovery Program. Li F. C. Sun Z. Hu D. K. Wang Z. W. (2013a). Crustal Structure and Deformation Associated with Seamount Subduction at the North Manila Trench Represented by Analog and Gravity Modeling. Mar. Geophys. Res. 34, 393–406. 10.1007/s11001-013-9193-5 Li Z. H. Xu Z. Q. Gerya T. Burg J. P. (2013b). Collision of Continental Corner from 3-d Numerical Modeling. Earth Planet. Sci. Lett. 380, 98–111. 10.1016/j.epsl.2013.08.034 Liao J. Gerya T. Malausà M. G. (2018). 3D Modeling of Crustal Shortening Influenced by Along-Strike Lithological Changes: Implications for Continental Collision in the Western and Central Alps. Tectonophysics 746, 425–438. 10.1016/j.tecto.2018.01.031 Liu S. Q. Zhao M. H. Sibuet J. -C. Qiu X. L. Wu J. Zhang J. Z. (2018). Geophysical Constraints on the Lithospheric Structure in the Northeastern South China Sea and its Implication for the South China Sea Geodynamics. Tectonophysics 742-743, 101–119. 10.1016/j.tecto.2018.06.002 Liu Z. Dai L. M. Li S. Z. Wang L. L. Xing H. L. Liu Y. J. (2021). When Plateau Meets Subduction Zone: A Review of Numerical Models. Earth-Science Rev. 215, 103556. 10.1016/j.earscirev.2021.103556 Magni V. Allen M. B. van Hunen J. Bouilhol P. (2017). Continental Underplating after Slab Break-Off. Earth Planet. Sci. Lett. 474, 59–67. 10.1016/j.epsl.2017.06.017 Manalo P. C. Dimalanta C. B. Faustino-Eslava D. V. Ramos N. T. Queaño K. L. Yumul G. P. Jr. (2015). Crustal Thickness Variation from a Continental to an Island Arc Terrane: Clues from the Gravity Signatures of the Central Philippines. J. Asian Earth Sci. 104, 205–214. 10.1016/j.jseaes.2014.08.031 Manea V. C. Manea M. Ferrari L. Orozco-Esquivel T. Valenzuela R. W. Husker A. (2017). A Review of the Geodynamic Evolution of Flat Subduction in Mexico, Peru, and Chile. Tectonophysics 695, 27–52. 10.1016/j.tecto.2016.11.037 Martinod J. Funiciello F. Faccenna C. Labanieh S. Regard V. (2005). Dynamical Effects of Subducting Ridges: Insights from 3-D Laboratory Models. Geophys. J. Int. 163, 1137–1150. 10.1111/j.1365-246x.2005.02797.x Mason W. G. Moresi L. Betts P. G. Miller M. S. (2010). Three- Dimensional Numerical Models of the Influence of a Buoyant Oceanic Plateau on Subduction Zones. Tectonophysics 483, 71–79. 10.1016/j.tecto.2009.08.021 McIntosh K. Lavier L. van Avendonk H. Lester R. Eakin D. Liu C. S. (2014). Crustal Structure and Inferred Rifting Processes in the Northeast South China Sea. Mar. Pet. Geol. 58, 612–626. 10.1016/j.marpetgeo.2014.03.012 McIntosh K. van Avendonk H. Lavier L. Lester R. Eakin D. Wu F. (2013). Inversion of a Hyper-Extended Rifted Margin in the Southern Central Range of Taiwan. Geology 41 (8), 871–874. 10.1130/G34402.1 Miller M. S. Gorbatov A. Kennett B. L. N. (2005). Heterogeneity within the Subducting Pacific Slab beneath the Izu-Bonin-Mariana Arc: Evidence from Tomography Using 3D Ray Tracing Inversion Techniques. Earth Planet. Sci. Lett. 235, 331–342. 10.1016/j.epsl.2005.04.007 Miller M. S. Kennett B. L. N. Lister G. S. (2004). Imaging Changes in Morphology, Geometry, and Physical Properties of the Subducting Pacific Plate along the Izu-Bonin-Mariana Arc. Earth Planet. Sci. Lett. 224, 363–370. 10.1016/j.epsl.2004.05.018 Morra G. Regenauer-Lieb K. Giardini D. (2006). Curvature of Oceanic Arcs. Geology 34 (10), 877–880. 10.1130/G22462.1 Ranalli G. (1995). Rheology of the Earth. London: Chapman & Hall. Rosebaum G. Giles D. Saxon M. Betts P. G. Weinberg R. F. Duboz C. (2005). Subduction of the Nazca Ridge and the Inca Plateau: Insights into the Formation of Ore Deposits in Peru. Earth Planet. Sci. Lett. 239, 18–32. 10.1016/j.epsl.2005.08.003 Ruh J. B. (2016). Submarine Landslides Caused by Seamounts Entering Accretionary Wedge Systems. Tera Nova 0, 1–8. 10.1111/ter.12204 Schmeling H. Babeyko A. Enns A. Faccenna C. Funiciello F. Gerya T. (2008). A Benchmark Comparison of Spontaneous Subduction Models – toward a Free Surface. Phys. Earth Planet Inter. 171, 198–223. 10.1016/j.pepi.2008.06.028 Seno T. Stein S. Gripp A. E. (1993). A Model for the Motion of the Philippine Sea Plate Consistent with NUVEL-1 and Geological Data. J. Geophys. Res. 98, 17941–17948. 10.1029/93JB00782 Sibuet J. -C. Yeh Y. C. Lee C. S. (2016). Geodynamics of the South China Sea. Tectonophysics 692, 98–119. 10.1016/j.tecto.2016.02.022 Simons W. J. F. Socquet A. Vigny C. Ambrosius B. A. C. Abu S. H. Promthong C. (2007). A Decade of GPS in Southeast Asia: Resolving Sundaland Motion and Boundaries. J. Geophys. Res. 112, B06420. 10.1029/2005JB003868 Sun Z. Stock J. Klaus A. the Expedition 367 Scientists (2018). Expedition 367 Preliminary Report: South China Sea Rifted Margin. Int. Ocean. Discov. Program. 10.14379/iodp.pr.367.2018 Sun Z. Zhong Z. H. Keep M. Zhou D. Cai D. S. Li X. S. (2009). 3D Analogue Modeling of the South China Sea: A Discussion on Breakup Pattern. J. Asian Earth Sci. 34, 544–556. 10.1016/j.jseaes.2008.09.002 Tang J. Chen L. Meng Q. (2020). The Effects of the Thermal State of Overriding Continental Plate on Subduction Dynamics: Two-Dimensional Thermal-mechanical Modeling. Sci. China Earth Sci. 63, 1519–1539. 10.1007/s11430-019-9624-1 Tao J. L. Dai L. M. Lou D. Li Z. H. Zhou S. H. Liu Z. (2020). Accretion of Oceanic Plateaus at Continental Margins: Numerical Modeling. Gondwana Res. 81, 390–402. 10.1016/j.gr.2019.11.015 Tetreault J. L. Buiter S. J. H. (2014). Future Accreted Terranes: A Compilation of Island Arcs, Oceanic Plateaus, Submarine Ridges, Seamounts, and Continental Fragments. Solid earth. 5, 1243–1275. 10.5194/se-5-1243-2014 Tetreault J. L. Buiter S. J. H. (2012). Geodynamic Models of Terrane Accretion: Testing the Fate of Island Arcs, Oceanic Plateaus, and Continental Fragments in Subduction Zones. J. Geophys. Res. 117, B08403. 10.1029/2012JB009316 Tu K. Flower M. F. J. Carlson R. W. Xie G. H. Chen C. Y. Zhang M. (1992). Magmatism in the South China Basin: 1. Isotopic and Trace-Element Evidence for an Endogenous Dupal Mantle Component. Chem. Geol. 97, 47–63. 10.1016/0009-2541(92)90135-R Turcotte D. L. Schubert G. (2002). Geodynamics. Cambridge, UK: Cambridge University Press. Ueda K. Gerya T. V. Sobolev S. V. (2008). Subduction Initiation by Thermal-chemical Plumes. PhysEarth Planet Inter 171, 296–312. 10.1016/j.pepi.2008.06.032 van Hunen J. van den Berg A. P. Vlaar N. J. (2002). On the Role of Subducting Oceanic Plateaus in the Development of Shallow Flat Subduction. Tectonophysics 352, 317–333. 10.1016/S0040-1951(02)00263-9 van Hunen J. van den Berg A. P. Vlaar N. J. (2004). Various Mechanisms to Induce Presentday Shallow Flat Subduction and Implications for the Younger Earth: a Numerical Parameter Study. Phys. Earth Planet. Interiors 146, 179–194. 10.1016/j.pepi.2003.07.027 Vogt K. Gerya T. (2014). From Oceanic Plateaus to Allochthonous Terranes: Numerical Modeling. Gondwana Res. 25, 494–508. 10.1016/j.gr.2012.11.002 Wang C. Y. Ding W. W. Schellart W. P. Li J. B. Dong C. Z. Fang Y. X. (2021). Effects of Multi-Seamount Subduction on Accretionary Wedge Deformation: Insights from Analogue Modelling. J. Geodyn. 145, 101842. 10.1016/j.jog.2021.101842 Wang L. L. Dai L. M. Gong W. Li S. Z. Jiang X. D. Foulger G. (2022). Subduction Initiation at the Solomon Back-Arc Basin: Contributions from Both Island Arc Rheological Strength and Oceanic Plateau Collision. Geophys. Res. Lett. 49, e2021GL097666. 10.1029/2021GL097666 Wang T. K. Chen M. K. Lee C. S. Xia K. Y. (2006). Seismic Imaging of the Transitional Crust across the Northeastern Margin of the South China Sea. Tectonophysics 412, 237–254. 10.1016/j.tecto.2005.10.039 Wang Y. J. Han X. Q. Luo Z. H. Qiu Z. Y. Ding W. W. Li J. B. (2009). Late Miocene Magmatism and Evolution of Zhenbei-Huangyan Seamount in the South China Sea: Evidence from Petrochemistry and Chronology. Acta Oceanol. Sin. 31, 93–102. (in Chinese with English abstract). 10.1109/CLEOE-EQEC.2009.5194697 Wu J. Suppe J. Lu R. Q. Kanda R. (2016). Philippine Sea and East Asian Plate Tectonics since 52 Ma Constrained by New Subducted Slab Reconstruction Methods. J. Geophys. Res. Solid Earth 121, 4670–4741. 10.1002/2016JB012923 Wu J. Suppe J. (2018). Proto-South China Sea Plate Tectonics Using Subducted Slab Constraints from Tomography. J. Earth Sci. 29 (6), 1304–1318. 10.1007/s12583-017-0813-x Yan Q. S. Shi X. F. Wang K. S. Bu W. R. Xiao L. (2008). Major Element, Trace Element, and Sr, Nd and Pb Isotope Studies of Cenozoic Basalts from the South China Sea. Sci. China Ser. D Earth Sci. 51, 550–566. 10.1007/s11430-008-0026-3 Yan Z. Y. Chen L. Xiong X. Wan B. Xu H. Z. (2021). Ocanic Plateau and Subduction Zone Jump: Two-Dimensional Thermo-Mechanical Modeling. J. Geophys. Res. Solid Earth 126, e2021JB021855. 10.1029/2021JB021855 Yan Z. Y. Chen L. Xiong X. Wang K. Xie R. X. Hsu H. (2020). Observations and Modeling of Flat Subduction and its Geological Effects. Sci. China Earth Sci. 63, 1069–1091. (in Chinese with English abstract). 10.1007/s11430-019-9575-2 Yang S. H. Li Z. H. Gerya T. Xu Z. Q. Shi Y. L. (2018). Dynamics of Terrane Accretion during Seaward Continental Drifting and Oceanic Subduction: Numerical Modeling and Implications for the Jurassic Crustal Growth of the Lhasa Terrane, Tibet. Tectonophysics 746, 212–228. 10.1016/j.tecto.2017.07.018 Yang T. F. Lee T. Chen C. H. Cheng S. N. Knittel U. Punongbayan R. S. (1996). A Double Island Arc between Taiwan and Luzon: Consequence of Ridge Subduction. Tectonophysics 258, 85–101. 10.1016/0040-1951(95)00180-8 Yeh Y. C. Hsu S. K. Doo W. B. Sibuet J. -C. Liu C. S. Lee C. S. (2012). Crustal Features of the Northeastern South China Sea: Insights from Seismic and Magnetic Interpretations. Mar. Geophys. Res. 33, 307–326. 10.1007/s11001-012-9154-4 Yeh Y. C. Sibuet J. -C. Hsu S. K. Liu C. S. (2010). Tectonic Evolution of the Northeastern South China Sea from Seismic Interpretation. J. Geophys. Res. Solid Earth 115, B06103. 10.1029/2009jb006354 Yoshida M. (2017). Trench Dynamics: Effects of Dynamically Migrating Trench on Subducting Slab Morphology and Characteristics of Subduction Zones Systems. Phys. Earth Planet Inter. 268, 35–53. 10.1016/j.pepi.2017.05.004 Zhang J. Li J. B. Ruan A. G. Ding W. W. Niu X. W. Wang W. (2020). Seismic Structure of a Postspreading Seamount Emplaced on the Fossil Spreading Center in Southwest Subbasin on the South China Sea. J. Geophys. Res. Solid Earth 125, e2020JB019827. 10.1029/2020JB019827 Zhao M. H. He E. Y. Sibuet J. -C. Sun L. T. Qiu X. L. Tan P. C. (2018). Postseafloor Spreading Volcanism in the Central East South China Sea and its Formation through an Extremely Thin Oceanic Crust. Geochem. Geophys. Geosyst. 19, 1–21. 10.1002/2017GC007034 Zhao Y. H. Ding W. W. Ren J. Y. Li J. B. Tong D. J. Zhang J. Y. (2021). Extension Discrepancy of the Hyper-Thinned Continental Crust in the Baiyun Rift, Northern Margin of the South China Sea. Tectonics 40, e2020TC006547. 10.1029/2020tc006547