Subduction Initiation at the Corner of Small Oceanic Basins

In Southeast Asia, emerging subduction zones often appear to begin at the corners of small oceanic basins, which have a triangular‐indenter continent–ocean boundary geometry. To investigate the influence of a triangular corner on subduction initiation, we performed a series of three‐dimensional numerical simulations with varying corner angles and base lengths. The results show that the apex of the corner constitutes the initial location of subduction, irrespective of the angle or the extent of the corner. Smaller angle corners are more likely to facilitate subduction initiation. At the same time, wide acute angle corners are difficult to form. Our findings suggest that triangular corner structures may facilitate subduction initiation in smaller basins; however, the role such corners in subduction initiation is limited in larger basins. Our results emphasize the importance of accounting for the three‐dimensional geometry of a subduction zone when examining its subduction dynamics and geological features.


Introduction
Plate tectonic activity is an essential prerequisite for Earth habitability and sets Earth apart from other planets (Korenaga, 2012).As the foundation of plate tectonics, understanding the initiation of subduction remains a significant challenge and frontier in the field of Earth science (Gurnis et al., 2004;Korenaga, 2013;Stern, 2004;Stern & Gerya, 2018;Toth & Gurnis, 1998).However, one of the most notable difficulties in subduction initiation research is the lack of "direct observations" (Mueller & Phillips, 1991).Regions that are actively undergoing subduction initiation, such as Puysegur, can very rarely be directly observed (Gurnis et al., 2004;Shuck et al., 2022).Therefore, the subduction initiation process is usually studied using various techniques such as numerical simulations, physical simulations, and petrology (Agard et al., 2007(Agard et al., , 2016;;Bercovici, 2003;Dong et al., 2022;Gerya et al., 2015;Ishizuka et al., 2011;Maunder et al., 2020;McKenzie, 1977;Niu et al., 2003;Zhou et al., 2020) because mature subduction zones have undergone long periods of tectonic evolution, resulting in the disappearance of most structural features dating from the initial stage of subduction.Therefore, it is important to observe and study subduction initiation in areas where subduction initiation is either in progress or has recently been completed (Hall, 2019;Shuck et al., 2022).
Southeast Asia is home to several small oceanic basins and subduction zones that are relatively young (Figure 1a; Hall, 1996Hall, , 2012;;Hall & Spakman, 2015;Lai et al., 2021;Li et al., 2023).Examples include the Cotabato and Sulawesi subduction zones, which were initiated during the Pliocene or Late Miocene (Kopp et al., 1999;Schlüter et al., 2001;Silver et al., 1983) and have subduction slabs with a maximum depth in the middle of the trench and a gradual decrease in length along both sides (Figures 1c and 1d).Other regions have ongoing subduction initiation, such as the Tolo Trough (Figure 1b), where a portion of the continental crust is migrating over the oceanic crust of the Northern Banda Sea (Hall, 2019;Husson et al., 2022;Silver et al., 1983).In areas such as the Sula Deep (Figure 1b), which has abnormally deep water and steep terrain, subduction initiation seems probable in the future (Hall, 2019).Despite differences in their subduction initiation stages, these regions share similar structural characteristics, with oceanic lithosphere subducting beneath either a relic arc or a continental fragment.Moreover, and more importantly, subduction initiation, either ongoing or imminent, is consistently located in the corners of these oceanic basins.
Subduction initiation is a complex process that is affected by multiple factors, including the density differences between the plates, thermal structure, and weak zones (Leng & Gurnis, 2015;Nikolaeva et al., 2010;Niu et al., 2003;Zhou et al., 2020).For example, Nikolaeva et al. (2010Nikolaeva et al. ( , 2011) ) conducted numerical simulations of subduction initiation in a passive continental margin and proposed that subduction initiation is facilitated by a thick crust and a thin, hot lithosphere.Similarly, Leng and Gurnis (2015) concluded that, based on numerical simulations, subduction initiation is easier in arc regions than in passive margins.The overriding plates in the small, mostly island arc, subduction zones in Southeast Asia are consistent with the results of these studies.However, these studies were based on two-dimensional (2D) numerical models.Some lateral structural geometry changes affecting subduction initiation may require three-dimensional (3D) numerical models.With the rapid  Seton et al. (2020).NBS indicates the North Banda Sea, and SBS indicates the South Banda Sea.(b) Bathymetric image of the North Banda Sea and surrounding regions, taken from Figure 7 in Hall (2019).(c) Subducted Celebes Sea Plate in the North Sulawesi and Cotabato trenches (the gray irregular polygons indicate the shape of the subducting plate projected onto the surface, while the red irregular polygons are reconstructions of the subducted plate).(d) Diagram of subducted plate reconstruction to the surface.The black line represents the subducted slab.The calculated length is shown on the surface as a red line.The slab dip along the C-D profile was not certain from earthquake events and was extracted from the Slab2 model (Hayes et al., 2018).development of computing power, more 3D models are now being used to study subduction initiation (Almeida et al., 2022;Gerya et al., 2015;Riel et al., 2023;Zhou et al., 2018Zhou et al., , 2020)).Marques et al. (2014) used three 3D models to investigate subduction initiation at a passive continental margin (where the angle is zero between the trace of the curved section of the margin and the z-axis of the model), as well as other models at different angles to the z-axis.Their results showed that a straight passive margin is more susceptible to subduction initiation than a curved passive margin.However, these models did not incorporate geometric shapes such as a "triangular corner," which refers to a specific location where the oceanic plate boundary geometry exhibits a change in direction, causing it to zigzag toward the continent and form an angle.Significantly, the small oceanic basins in Southeast Asia differ significantly from the long, generally straight margins of the Atlantic.In Southeast Asia, the continent-ocean boundaries are typically much narrower, formed more quickly, and are much steeper (Hall, 2019).Currently, there are no numerical models available to validate the observed subduction-prone corners in Southeast Asia.

Model Setup
We employed the massively parallel finite element code ASPECT (Bangerth et al., 2021;Heister et al., 2017;Kronbichler et al., 2012) to construct and execute our 3D numerical geodynamic model (Figure 2a).The left side of the model represents the oceanic plate, composed of the oceanic crust (10-km thick) and the suboceanic lithospheric mantle.The right side of the model represents the overriding plate, comprising the upper crust (25-km thick), lower crust (15-km thick), and lithospheric mantle (40-km thick).This configuration creates a substantial density contrast between the overriding plate and the oceanic plate (Figure 2b).We adopted this specific setup for the overriding plate model based on previous studies (Marques et al., 2014;Nikolaeva et al., 2010), allowing the model to achieve completion within a short time.Detailed information concerning the initial model design, equations, and parameters can be found in Text S1 in Supporting Information S1.
A 10-km-thick weak zone between the oceanic plate and the overriding plate is included in our model to ensure that subduction initiation can be completed.Typically, the boundary between the oceanic plate and the passive continental margin is not vertical in the z-direction and part of the oceanic plate extends below the continental plate (Marques et al., 2013(Marques et al., , 2014;;Nikolaeva et al., 2010Nikolaeva et al., , 2011)).Nevertheless, for the small basins in Southeast Asia, the overriding plates may closely resemble the lithosphere of relic arcs.Therefore, we adopted a similar approach to that of Leng and Gurnis (2015) by vertically dividing the plates according to their weak zones; this allowed us to isolate just the impact of the corner geometry on subduction initiation.
In our models, the geometry of the corner is defined by the angle α and width D of the base extent of the corner (Figure 2c).Model 0 (a straight ocean-continent boundary, invariant in the z-direction) serves as the reference model to evaluate the influence of the corner.Models 1, 2, and 3 have a width D of 200 km, with α values of 120°, 90°, and 60°, respectively.Models 4, 5, and 6 have a width D of 115.4 km, with α values of 120°, 90°, and 60°, respectively.Model 7 shares the same apex positions as Models 1 and 5, with an α value of 60°.Models 8, 9, and 10 have a width D of 346.4 km, with α values of 120°, 90°, and 60°, respectively.

Results
In this study, we designated the time at which the sinking oceanic slab reaches a depth of 200 km as the completion time of subduction initiation (Figure 3a; Figure S1 in Supporting Information S1), serving as a metric to evaluate the level of difficulty for each model to achieve subduction initiation.Because the temporal resolution of our model output is 0.1-Myr intervals, there may be an associated error of ±0.1 Myr in the recorded time (Table 1).
In comparison to the reference model (Model 0), characterized by a straight ocean-continent boundary, the majority of the models featuring a corner structure demonstrated accelerated subduction initiation.Only Model 1 exhibited a longer initiation period than the reference model.Note that, because of the marginal error of 0.1 Myr, the subduction initiation times in Models 4 and 8 can be approximated as being similar to those of the reference model.Interestingly, Models 1, 4, and 8 share the same angle of 120°.
In the simulations involving Models 1, 2, and 3, all featuring a corner structure with a consistent base extent of 200 km, we observed that a larger angle (α) was correlated with extended durations for completing subduction initiation.This trend was also evident in other models with varying D values, such as Models 4, 5, and 6, as well as Models 8, 9, and 10.Accordingly, we conclude that, for corner structures sharing the same base extent, a smaller angle tends to lead to quicker subduction initiation.
Despite these results, no consistent pattern emerged for models with the same angle but differing base lengths.For example, when considering an angle of 120°, Model 1, with D = 200 km, required 5.0 Myr to achieve initial subduction.In contrast, Model 4, with a smaller corner, and Model 8, with a wider corner, both exhibited shorter subduction initiation times.
The results for models featuring a 90°angle diverged from those with a 120°angle.Notably, Model 2, with a corner base extent of 200 km, exhibited a higher probability of subduction initiation than both the narrower corner (Model 5) and the wider corner (Model 9), which required more time to complete initial subduction than Model 2. Nevertheless, all models with a 90°angle accomplished initial subduction more rapidly than Model 0.
In the case of models with an angle of 60°, analogous to those with an angle of 90°, the model (Model 3) with D = 200 km completed subduction initiation faster, whereas both of the narrower corners (Models 6 and 7) and the wider corner (Model 10) exhibited lengthier subduction initiation times.However, the time difference was not as pronounced as observed in models with a 90°angle.Considering the potential error of 0.1 Myr, Models 3 and 10 are reasonably similar.Consequently, the relationship between the time required for initial subduction completion and the base extent D is not linearly dependent given a constant corner angle.
Analyzing models sharing the same apex position (h), yet varying in angle, we reached similar conclusions.For example, Models 1, 5, and 7, possessing identical apex positions, demonstrated that a smaller angle leads to a briefer duration for initial subduction.However, within the context of models with the same angle, the discrepancy in the subduction initiation time arising from different h values did not adhere to a linear pattern.

Influence of the 3D Corner Geometry on Subduction Initiation
In a 3D model, the magnitude of the vertical downward force (the y-component) can be amplified or diminished given the presence of a horizontal component (the z-component).This component, a representation of the third spatial dimension lacking in a 2D model, influences the distribution of forces within the 3D framework.
For the first few timesteps (less than 1 Myr), the continental plate undergoes rapid vertical movement to reach isostatic equilibrium (Figure 4b; Figures S2b and S3b in Supporting Information S1).Subsequently, the primary driving force, originating from the density contrast between the oceanic and continental plates (Figure 2b), causes the overriding plate to overthrust onto the oceanic plate, inducing downward bending of the latter.Particularly at the apex of the corner, a focusing effect facilitates the dipping of the oceanic plate.As subduction initiation progresses, the leading edge increases its buoyancy contrast with the surrounding mantle, intensifying the bending and the downward pull of the plate.As subduction initiation progresses, the corner area on the surface gradually decreases (Figure 4b; Figures S2b  and S3b in Supporting Information S1).However, the extent of this reduction varies among the models with different angles.Models with smaller initial angles experience greater reductions in the extent of the corner (D) and in the reduction of h (∆h).In models with large angle corners, the force component in the z-direction is smaller than that in the x-direction.While the apex position of the corner has a force in the z-direction on both sides that pushes the oceanic plate downward, this force is too small to change the deformation of the overriding plate in the corner.As a result, the change in the length h is relatively small and the overriding plate moves toward the oceanic plate along the x-direction, with the continent-ocean boundary shape remaining largely unchanged (Figure 4b).Therefore, large-angle corners have little effect on subduction initiation because they lose the focusing effect of the corner as nearly the entire margin geometry becomes triangular.According to our numerical simulations, the corner does not promote subduction initiation when the angle of the corner exceeds 120°.
In the case of a small angle corner, the component of the force in the z-direction is sufficiently large to induce deformation of the overriding plate.Consequently, the oceanic lithosphere at the apex of the corner is compressed by the overriding plate on both sides, causing the apex position of the corner to gradually close.This leads to the oceanic plate at the apex sinking rapidly into the asthenosphere.As the plate reaches asthenospheric depths, the buoyancy contrast and viscosity contrast between the subducting plate and the surrounding mantle increase.Consequently, slab pull becomes increasingly significant.Therefore, the sooner the oceanic plate at the corner enters the asthenosphere, the earlier the entire subduction system will form (Figure 3a; Figure S1 in Supporting Information S1).

Initial Subduction Starts at the Corner Point
In all models, the apex position of the corner is the site at which the model plate first begins to bend downward, irrespective of the angle or extent of the corner.The subducting slab reaches a depth of 200 km first at the midsection of the plate (as shown in Figure 3b).This is consistent with the conditions observed in the Celebes Sea, where the subducted slab has reached the deepest point at middle of the trench both in the Cotabato and North Sulawesi subduction zones.
The current plan views of the subduction zones appear to be relatively flat without prominent corners (gray areas in Figure 1c).However, this does not indicate the shape of the plate prior to subduction.On the basis of the Slab2.0model, we can calculate the length of the subducted slab along the trench-normal direction.Then, we can reconstructed the slab to its pre-subduction shape (Figure 1d).We can see that there is a triangle corner at the front end of the reconstructed plate for both Cotabato and North Sulawesi (red areas in Figure 1c).The longest slab, corresponding to the location where the subduction initiation first starts, is the apex of the corner.This is consistent with our simulation results, which show that subduction initiation starts at the apex of the corner.
Our modeling results also correspond to the conditions observed in the North Banda Sea.The region's deepest points correspond to corner apex positions (Figure 1b).The northeastern corner of this region (the Sula Deep) exhibits depths exceeding 5,800 m, characterized by slope gradients of up to 13° (Hall, 2019).These gradients are likely related to slab pull at depth resulting from roll back in the Banda Sea (Husson et al., 2022).The Tolo Trough, located at the western edge, exhibits average slope gradients of 7-8°near its deepest segment.Rudyawan and Hall (2012) have documented thrusting in this region, consistent with downslope movement of a thrust sheet and the displacement of continental crust onto the North Banda Sea oceanic crust.
Both numerical simulations and geological observations consistently indicate that subduction initiation predominantly takes place at the corner apex, irrespective of the corner angle or extent.This phenomenon arises from the inherent property of the corner geometry, which concentrates stress at its vertices, particularly during the initial stage of subduction (Figure 4a; Figures S2a, and S3a in Supporting Information S1).The apexes of corners, where stress concentrates, are more susceptible to deformation than other areas.In addition, the stress concentration in the corner of a corner with a smaller angle is greater; this is the primary reason why the corner model with a smaller angle is more likely to speed up subduction initiation.
Although an acute angle corner geometry can easily trigger subduction initiation, this structure is unstable because it is easily compressed by both sides of the plate (Figure 4).Consequently, forming a wide range of such sharp corners in actual geological cases is challenging.If these structures were to form during oceanic plate development, before reaching significant size, intense stresses would rapidly compress the corner area, particularly at the apex, potentially leading to subduction.As a result, such large acute angle corner structures are not observed on the Earth's surface at present.
The promotion of subduction initiation by a corner primarily relies on the negative buoyancy of the plate bending downward at the apex position.In smaller basins (hundreds of kilometers in size), inverting the short margin necessitates a smaller force, allowing a small acute or right-angled corner to facilitate subduction initiation.However, in larger basins (several thousands of kilometers in scale), inverting the entire long margin demands a larger force, indicating the need for a large acute angular corner.Unfortunately, such a large corner formation is improbable.Therefore, while the triangular corner structure may facilitate subduction initiation in smaller basins, its role in larger basins is limited.

Model Limitations
In this study, all of the simulations accounted for the weak zone at the ocean-continent boundary, with different corner geometries corresponding to various weak zone configurations.However, these models are idealized.In geological settings, the strength of the ocean-continent boundary relies on inherited structures and the degree of serpentinization, which may not uniformly exhibit low strength across the entire boundary.Consequently, subduction initiation is more likely to occur where weak strength areas are present (Zhou et al., 2020).
While our models were designed to simulate spontaneous subduction, it is crucial to acknowledge that small basins, often situated amidst numerous active plates, can experience influences from surrounding plate activities during subduction initiation (Dong et al., 2022).Quantifying external forces acting in different directions within a 3D model is complex.Therefore, we intentionally omitted external forcing, focusing solely on assessing the impact of the corner geometry on subduction.Further research involving more complex structures and boundary conditions closer to those of natural subduction zones is necessary to enable a comprehensive understanding of subduction initiation dynamics.

Conclusions
In this study, we simulated a series of numerical models to investigate the relationship between the triangular corner geometry and subduction initiation.The results of the study suggest the following conclusions.
1.In small oceanic basins, subduction initiates at the apex of the corner, regardless of its size or angle.2. A small angle corner in an oceanic basin is favorable to subduction initiation.
3. An corner with an angle of more than 120°does not facilitate subduction initiation and may instead hinder it.

Figure 1 .
Figure 1.(a) Oceanic basins and subduction zones in Southeast Asia.The ages of the oceanic basins are taken from Seton et al. (2020).NBS indicates the North Banda Sea, and SBS indicates the South Banda Sea.(b) Bathymetric image of the North Banda Sea and surrounding regions, taken from Figure 7 in Hall (2019).(c) Subducted Celebes Sea Plate in the North Sulawesi and Cotabato trenches (the gray irregular polygons indicate the shape of the subducting plate projected onto the surface, while the red irregular polygons are reconstructions of the subducted plate).(d) Diagram of subducted plate reconstruction to the surface.The black line represents the subducted slab.The calculated length is shown on the surface as a red line.The slab dip along the C-D profile was not certain from earthquake events and was extracted from the Slab2 model(Hayes et al., 2018).

Figure 2 .
Figure 2. (a) Initial stage showing the triangular corner structure to simulate curved margins, the base length of the corner (D), and the angle of the corner (α).A zoom-in of the computational mesh shows that the minimum resolution is ∼3 km.(b) Plate strength and density contrast between the oceanic and overriding plates.(c) Corner geometries in different model experiments.

Figure 3 .
Figure 3. (a) Final results of Models 0-3 showing only the oceanic crust and suboceanic lithospheric mantle.Models 1-3 have the same D (200 km) and different angles (α).(b) Sections through end results (at 3.2 Myr) for Model 2 at y = 0 km (front), y = 200 km (middle), and y = 400 km (back).The oceanic slab reaches a depth of 200 km faster in the middle than in the front and back.The colors indicate the composition, as defined in Figure 2a.

Figure 4 .
Figure 4.Sections through models at z = 0 km.(a) Shear stresses (τ xy ) on the surfaces of Models 1, 2, and 3 at the initial stage.(b) Evolution of the corner geometry.The dashed lines represent the lengths of D and h in the initial stage.∆h represents the decrease in the length of h.The black arrows indicate the velocity.

Table 1
Description of Numerical Experiments With Varying Length and Angle of the Corner