Introduction

Oceanic subduction is one of the most prominent tectonic processes on the Earth’s surface. Theoretically, from ~20 to 40 Myr after formation, the oceanic lithosphere has cooled and been dense enough to be gravitationally unstable and thereby subductable1,2. Thus, the wide distribution of oceanic subduction on Earth appears normal. However, >50% of the present-day Atlantic Ocean floor is older than 40 Myr without widespread subduction3, and some oceanic blocks, like in the eastern Mediterranean, have remained at Earth’s surface for >300 Myr4. According to previous studies, the tendency toward subduction of aged oceans can be suppressed by the bending resistance itself and, more importantly, by the frictional resistance of adjacent plates (i.e., a potential upper plate)5,6,7,8. Therefore, before subduction zone initiation (SZI)9, some precondition(s) are required to break the quasi-balanced state between the future upper and lower plates.

Based on previous studies, before SZI, the coupling between the oceanic plate and its surroundings should be broken first by the development of a mechanically weak zone in between, e.g., due to intra-oceanic transform faulting5,10,11, new oceanic ridge spreading10, passive margin rifting12, or the passive margin’s lithospheric delamination13. Then, the appearance of vertical or horizontal-dominated external forces is found necessary to further trigger SZI (Fig. 1), e.g., due to large-scale mantle convection14,15,16,17,18, tectonic loading at the passive margin2,19, or trenchward contraction (or plate push) due to surrounding plate interactions6,8,20,21,22,23,24. Systematically, oceanic subduction triggered by vertical or horizontal-dominated forces is expected to show differences in the time sequence between the initial burial of oceanic materials and the rifting of the upper plate—nearly contemporary (Fig. 1a, b) versus burial before rifting (Fig. 1c, d)8,25,26. Therefore, the early-stage evolution of subduction appears key to understanding the triggering mechanisms for modern oceanic subduction.

Fig. 1: Two main forcing scenarios for SZI development.
figure 1

a, b Vertically forced scenario: After the thinning and weakening of the future upper plate, the oceanic plate’s negative buoyance leads to SZI10. Note that tectonic loading may also aid SZI, and like negative buoyancy, its direction is mainly vertical2,19. c, d Horizontally forced scenario: Horizontal force (or plate push) overcomes the resistance from the future upper plate and triggers SZI8,25,29. e Stage of mature subduction. This stage features horizontal mantle wedge flow, sparse arc magmatism, and serpentinites developed in the subducting lithospheric mantle or eroded from the forearc mantle as the slab’s main dehydration components29,37,39,47. The decoupling zone could be the junction between compositionally/rheologically different tectonic blocks that favor strain localization8,66. The initial condition in this cartoon figure specifically targets an Izu–Bonin–Mariana-type SZI development we mainly simulated in this study. The reader should check more SZI scenarios in the literature, e.g., Crameri et al.9.

Full size image

Due to the prolonged subduction erosion or entrainment27,28, in most subduction zones, the geological records for the initial subduction evolution have been largely destroyed and challenging to discriminate8. Fortunately, thanks to the International Ocean Discovery Program (IODP) expedition drilling, as well as diving and dredging, a complete SZI magmatic stratigraphy has been recognized along the Cenozoic Izu–Bonin-Mariana (IBM) forearc (Fig. 2)29,30,31. Here, SZI occurred at ~52 Ma31, roughly when the Pacific plate-mantle system was reorganized32,33, and the Indo-Eurasia continental collision started to get intensive (Fig. 2a)34,35,36. Specifically, during the first ~7–8 Myr after SZI, the contribution of slab materials to subduction-related magmatic sources varied from nearly zero in earliest forearc basalts (FAB; 52–51 Ma), to primarily lower-ocean crust contributions in the subsequent low silica boninites (51–50 Ma), to sediment and upper-ocean crust contributions in later-stage high silica boninites (50–46 Ma; Fig. 2c), followed finally by growing contributions from the lithospheric mantle of the subducted slab (i.e., fluid/melt generated via serpentinite decomposition) in mature-arc magmas (<44 Ma)29,30,37,38. Therefore, magma geochemistry in the IBM subduction zone could provide the necessary fingerprints for constraining the initial scenarios for oceanic subduction29,30 (Figs. 1 and 2). In addition, studies of this area provide some vital information about space-time evolution during SZI: (1) the newly formed seafloor due to syn-SZI rifting is ~300 km wide39; (2) the maximum trench-away migration of the magmatic center in the upper plate after the SZI is >80 km (Fig. 2b); and (3) the evolution from the SZI to a mature subduction zone took ~7–8 Myr (Fig. 1e)39,40,41. These observations provide primary constraints for testing forward models of modern oceanic subduction.

Fig. 2: Tectonic background and key observations of the Izu–Bonin-Mariana (IBM)-type subduction zone.
figure 2

a Tectonic background for SZI in the Izu–Bonin-Mariana-type subduction system3,32,34. At ~52 Ma, the Pacific plate is suggested to be separated from the Philippine plate by a weak transform fault8,9,14,62,63, the motion direction of the Pacific plate rapidly changed32,33, and the continental plates with gray colors were roughly drifting towards the IBM, i.e., Australia, India, Arabia, and the terranes in present-day Southeast Asia. The brown region marks one possible reconstruction of the subducted block(s) between India and Eurasia34. b Map of the IBM sample localities in (c) (after Ishizuka et al.30). The subduction evolution is nearly synchronous throughout the >2500 km long trench (see inset)31. c Pb isotopes of the IBM igneous rocks reveal changing subduction inputs with time30,46. △8/4 (Pb) is the displacement from the Northern Hemisphere Reference Line (NHRL) in 208Pb/204Pb–206Pb/204Pb space67. The forearc basalts’ mantle source is similar to the Philippine Sea Mid-Ocean Ridge Basalts (MORB), without subduction inputs. Later Expedition 352 low silica boninites bear clear subduction inputs from slightly altered Pacific Oceanic crust, probably the lower crust29,46. Subsequent Expedition 352, Chichijima and Mukojima high silica boninites contain a substantial contribution from the marine sediments30,46,68. Marker meanings are the same as in (b).

Full size image

As most previous studies focused on the triggering mechanisms for SZI, a process that typically lacks observational constraints25, a comprehensive interpretation of the large-time-space observations above (Fig. 2) seems necessary to clarify the general evolution of modern oceanic subduction (Fig. 1). Therefore, we utilized the method of forward geodynamic and geochemical modeling, through connecting the time-space-discrete magmatic records (Fig. 2), to search for a best-fit scenario where IBM-type subduction tends to develop.

Results

Forward modeling

Following the initial conditions in previous studies11,42, we designed two types of numerical models (with vs. without horizontal force, i.e., trenchward contraction Fig. 3, corresponding to the two scenarios in Fig. 1) to explore the infancy-to-maturation evolution of modern oceanic subduction, especially targeting the IBM-type subduction zones where the geological records are similar along the entire length of the >2500 km long forearc (Fig. 2b)31. The initial and boundary conditions of numerical models are in Fig. 3 (model parameters are in Supplementary Table S1). When modeling the behavior of multi-phase melting, we treat the effects of melt-extraction-linked depletion and source petrology (crust and asthenospheric or lithospheric mantle) on the solidi (Supplementary Fig. S1). We use the phase diagrams for marine sediments, oceanic crust, and depleted mantle from Kimura (2017) to calculate the dehydration behaviors of slab materials (Fig. 4)34,43. We assume that the melting of subducted sediments and crust is crudely synchronous with their dehydration, as the presence of a free water phase should effectively lower their solidi and enhance melting44. In the numerical models, we define the moment of SZI as when the tip of the oceanic slab is deeper than 180 km deep, after which the subduction process becomes self-driven. The online method section provides a specific description of the tested variable ranges and a summary of the general features of the two model types. Our Supplemental Information provides more details as to the modeling method.

Fig. 3: Model setup.
figure 3

Initial and boundary conditions for the first-type (a) and second-type model (b) in Fig. 1. The insects in (a) demonstrate the material layering of the future upper and lower plates, while those in b present the initial rheological structures of both domains. Lith. Mantle- Lithospheric Mantle. In all numerical models, the 2D numerical box is 2500 km wide × 600 km deep42. The initial oceanic plate is 1800-km long. The finest resolution is ~5 km, where the mesh node is within 50 km from the nearest material interface, and it gradually increases to ~50 km, where the node is >190 km away from the closest material interface. Following previous studies, we assume that a preexisting weak zone (light green) decouples the further lower and upper plates8,11,14,42. The surface and bottom temperatures remain at 0 and 1350 °C, respectively.

Full size image
Fig. 4: Water contents of slab materials and the slab surface temperature in representative models.
figure 4

a–c Water contents of slab of lithospheric mantle (a), slab crust (b), and sediments (c) at different temperature-pressure conditions43. Lines with age demonstrate the evolution of slab surface temperature in the first-type model (Run 10). d–f The plotting habits are the same as in (a–c), except that slab surface temperature is from the second-type model (Run 71).

Full size image

A best-fit scenario of the IBM-type subduction

As mentioned earlier, the time-space pattern of magmatic distribution and geochemistry can be key to understanding the evolution of modern oceanic subduction45 (Fig. 2). By comparing the melting behavior in the forearc region and the dominant dehydration material(s) in the subducted slab with observations (Fig. 2), we find that the second-type model (Fig. 1c–e) with trenchward contraction can more effectively explain IBM-type subduction evolution (Fig. 5).

Fig. 5: Representative numerical model of the second-type scenario in Fig. 2 (Run 71, Fig. 9).
figure 5

a–h Model evolution regarding material, temperature (gray isolines), and velocity fields (gray arrows). Blue arrows illustrate the syn-SZI rifting locus in the upper plate, which also migrates towards the trench but slower than the melting locus. Material type: light blue—the ocean’s lithosphere mantle, green- lower oceanic crust, orange- upper oceanic crust (or altered oceanic crust), purple—sediments. i Volume evolution of scraped slab materials (demonstrated by the relative frequency of different material tracer types). Red and blue arrows denote the scraped materials’ sudden volume increase or decrease. The inset demonstrates the volume evolution of scraped slab materials in 30 Myr. j Time-space distribution of asthenospheric mantle melts. The red arrow demonstrates the migrating direction of mantle melts. k Dehydration behaviors of slab materials34,43. The modeled melts/fluids are shown for every 10 km along the horizontal direction, with the symbol size representing the number of tracers. A complete evolution of this model is in Supplementary Movies 1 and 6.

Full size image

Specifically, due to the trenchward contraction, the slab materials (i.e., marine sediments and upper crust) are scraped off after SZI, forming a small accretionary wedge in the forearc region (Fig. 5b, i). Therefore, the materials that dehydrate, which contribute to syn-SZI magmatism in the first ~2 Myr after SZI (~1.2 Myr), are mainly the lower (or slightly altered) slab crust, i.e., inputs for the Expedition 352 low silica boninites (Fig. 5k vs. 2c).

After SZI, the slab rolls back quickly (Fig. 5c), which causes upper-plate rifting and migration of the rifting locus away from the trench (Fig. 5d–g). Meanwhile, the net distance between the locus of mantle melting and the trench also increases (Fig. 5d–g, j). As the rifted and thinned upper plate weakens, the initially accreted slab materials can be entrained by subduction (Fig. 5c, i) and can dehydrate and contribute to the melting products of slab crust, i.e., inputs for high silica boninite (Fig. 5k vs. 2c).

As the rollback rate reduces and the rifted lithosphere solidifies, the rifting event ceases at ~8 Myr after SZI (Fig. 5g and Supplementary Movie 6). At this time, the center of mantle melting migrates ~100 km (Fig. 5j)30. Later, the convective flow’s direction in the mantle wedge becomes mainly horizontal (Fig. 5g), except for being occasionally disrupted by small-scale mantle drips (Fig. 5h). Because, before this time, a mature-arc-like mantle wedge flow does not exist, the formation of the earliest magmas (52–48 Ma FAB and boninites) leaves a localized mantle residuum whose depletion increases with melt extraction29,46. In contrast, after active corner flow develops in the mantle wedge region (Fig. 5g), the mantle source of the later magmas (~ 44 Ma Hahajima early arc magma) becomes more fertile, deviating from the boninites’ residual source mantle30. Subsequently, due to the solidification of the rifted lithosphere, the slab’s surficial materials are scraped off again at ~10 Myr after SZI (Fig. 5h, i), which later are consumed by the following subduction erosion (see inset in Fig. 5i). Together with the secular cooling of the slab surface, the primary materials dehydrating on the slab evolves to become hot lithospheric mantle (Fig. 5k). Therefore, the magmas of this stage (the Hahajima lavas) mainly show the Pb isotope character of unaltered oceanic crust30, close to those of the modern Izu arc volcanism37,47. At this moment, the model enters the mature stage of oceanic subduction40.

As a comparison, the average subduction rate of the initial oceanic plate \({{{{{{\rm{\upsilon }}}}}}}_{{{{{{\rm{ave}}}}}}}\) of a first-type model, Run 10, is like that of Run 71 (Fig. 6). However, after SZI, because the >300-km long tip portion of the initial oceanic plate suddenly sinks into the deep mantle, a comparably broad mantle window opens, where the upwelling asthenospheric mantle melts rapidly to create new lithosphere (Fig. 7a–j). Meanwhile, mantle magmatism and slab dehydration happen only before the solidification of this region (within ~3 Myr after the SZI) (Fig. 7l–m). During this brief stage, the magmatic locus migrates >300 km, unlike the IBM observation (Fig. 7l)30. In addition, since slab burial happens nearly synchronously with rifting, at the beginning of the SZI episode, all slab materials are observed to dehydrate (Fig. 7k), contrasting with our other observations (Fig. 2c).

Fig. 6: Sensitivity tests for the first-type model (Fig. 1).
figure 6

a The time needed for the subduction zone initiation (SZI) to occur. The SZI is identified when the left-most slab surface (i.e., the reference point in Fig. 3a) is deeper than the upper plate’s bottom. b The time when slab tearing happens (e.g., Fig. 10). nan demonstrates that slab tearing does not happen in the test. c The total time needed for the subduction of the initial oceanic plate (Fig. 3a). Because subduction re-initiation does not easily happen in this scenario (Supplementary Movie 2), nan denotes that the initial oceanic plate’s subduction does not finish within 100 Myr. d Average subduction rate. The average subduction rate is calculated only in models where the initial oceanic plate completely subducts. Rx is the xth model, and the number after “−“ denotes the corresponding time. \({{{{{{\rm{\rho }}}}}}}_{{{{{{\rm{AM}}}}}}}\) is the compositional density of the asthenospheric mantle. \({{{{{{\rm{Age}}}}}}}_{{{{{{\rm{slab}}}}}}}\) is the thermal age of the initial oceanic plate69. \({{{{{{\rm{\sigma }}}}}}}_{{{{{{\rm{slab}}}}}}}\) is the maximum yielding stress for the subducted slab (also see Supplementary Table S1). \({{{{{{\rm{\eta }}}}}}}_{\min }\) is the minimum viscosity (in Pa∙s) for the region behind the initial oceanic plate, e.g., a mid-oceanic ridge or newly formed oceanic lithosphere. In this figure, \({{{{{{\rm{\eta }}}}}}}_{\min }\) is in the form of lg (\({{{{{{\rm{\eta }}}}}}}_{\min }\)).

Full size image
Fig. 7: Representative numerical model of the second-type scenario in Fig. 2 (Run 10, Fig. 6).
figure 7

a–j Model evolution regarding material, temperature, and velocity fields (Initial condition in Fig. 3). k Volume evolution of scraped slab materials. l Time-space distribution of asthenospheric mantle melts. m Dehydration behaviors of slab materials. Plotting habits are the same as in Fig. 5. A complete evolution of this model is in Supplementary Movie 5.

Full size image

Discussion

Following the practices in previous studies14,42, we set a free-slip boundary condition at 600 km deep (Fig. 3). With this type of model setup, the simulated \({{{{{{\rm{\upsilon }}}}}}}_{{{{{{\rm{ave}}}}}}}\) cannot strictly mimic a natural condition where the lower mantle is not strong enough to prohibit slab sinking permanently. However, given that most subducted slabs in the IBM subduction zone are found to flatten upon the lower mantle to varying degrees48, and the total modeling time here is mostly ≤30 Myr, e.g., a potential duration for slab stagnation49—this model setup may still crudely fit the specific research area (Fig. 2). In addition, as the slab’s surface temperature in the second-type model gradually evolves to the observed ones (Fig. 8a), and the stress pattern within the slab is like that of most IBM slab portions (Fig. 8b vs. 8d), our numerical models appear consistent with the first-order observation of modern oceanic subduction.

Fig. 8: Comparison between model results and seismology-based studies.
figure 8

a Evolution of slab surface temperature. The shaded area demonstrates the slab surface temperature from global subduction zones70. Dashed lines are from Run 10 (the first scenario in Fig. 1). Solid lines are from Run 71 (the second scenario in Fig. 1). b Stress state within the IBM slabs71,72. c, d Vertical deviatoric stress of the subduction zone from Run 10 and Run 71, respectively, when slab stagnation happens.

Full size image

Within the scope here, the absence of shallow slab materials’ contribution (i.e., sediments and upper/altered oceanic crust) in the magma source during the first ~2 Myr after SZI29,46 is mainly due to trenchward contraction that scrapes these materials away (Fig. 5b, i); while the involvement of these materials later in the process is because that the weak rifted upper plate can no longer prevent their entrainment into subduction (Fig. 5c, i). The absence of sediment contributions in the earliest magmas with subduction inputs, i.e., low silica boninites, could also be due to their buoyancy or to upwelled mantle materials preventing their entrainment11,14. However, in these scenarios, the upper slab crust should still melt no later than the lower slab crust, which is inconsistent with previous studies29. In addition, the edge of the oceanic basin could also have been rejuvenated due to hot-spot-related upwelling, in which case its surface could rise to avoid receiving sediments, and its altered crust could be covered by fresh magma50. However, in this case, the slab density can also be reduced to lower the buoyancy contrast between the ocean and its adjacent plates51, and thus, the external force would be increasingly important to trigger SZI7; besides, due to the sparse distribution of reported hot spots, this scenario cannot easily explain the lack of sediment signal throughout the >2500 km long IBM trench46.

In the case of the IBM subduction, the trenchward contraction could result from the closure of the Tethys Ocean22, during which episodic terrane accretion and continental collision34,35 seem to have intensively enhanced plate interactions in the western Pacific region36,52 (Fig. 2a). Further, mantle suction, due to the rapid mantle downwelling following the previous avalanche of paleo-Pacific slabs could also contribute to the initial contraction15,16,17. As previous studies of plate reconstruction and global modeling reveal, at the time of SZI in the IBM system, the underlying mantle flow was most likely towards the Asian margin3,15. Therefore, if this external force had dominated, its effect would also be like the trenchward contraction assumed in the simplified models of this study (Fig. 3). We emphasize that the model-implied trenchward contraction could be sustaining or only transitory (~1–2 Myr) to break the state of quasi-balance between the future upper and lower plates (Figs. 3 and 9 and Supplementary Fig. S2). After that, the slab’s negative buoyancy will drive the subduction process and its active magmatism10,53, and sometimes, the SZI may even transfer to a vertically-driven stage20.

Fig. 9: Sensitivity tests for the second-type model (Fig. 1).
figure 9

a The time needed for the subduction zone initiation (SZI) to occur. The SZI is identified when the left-most slab surface (i.e., the reference point in Fig. 3b) is deeper than the upper plate’s bottom. “Always push” marks the runs where the added horizontal velocity for the initial oceanic plate is active throughout the model, while “2 Myr Push” demonstrates that it is only active during the first 2 Myr. b The time when slab tearing happens. c The total time needed for the subduction of the initial oceanic plate (Fig. 3a). d Average subduction rate. Plotting habits and abbreviations are the same as in Fig. 6.

Full size image

Conclusion

The modeling attempt here for the initial ~10–30 Myr evolution of deep subduction reveals that (1) during the earliest, non-magmatic stage of IBM subduction, trenchward contraction should have been active, which explains the absence of surficial slab material contributions in the earliest boninites; (2) the syn-SZI rifting event persists for ~8 Myr, during which the rifting and magmatic locus migrates ~100–150 km, and the rifted region is can be up to 300 km wide; and (3) subduction enters the mature stage after the rift’s solidification, at ~10 Myr after SZI—all features that mimic observations in the subduction evolution of the IBM system (Fig. 2 vs. 5). In addition to IBM, a similar evolutionary pattern has also been reported in the Cretaceous Oman subduction zone, another representative research area of modern deep oceanic subduction25,26,54. Therefore, the SZI development revealed in this study may also be universally instructive.

Methods

We use a code written in MATLAB (LL_AM_TMM) to conduct the 2D numerical experiments55,56,57,58,59. The code is based on the Lagrangian-type finite element code MILAMIN60. It includes a free surface on the top (a surface with no shear or normal stress)61 (Fig. 3) that has been benchmarked with a free subduction study55,57 (see the process of normal faulting in the bent oceanic plate in Supplementary Movie 1). The triangular mesh is adaptively generated/regenerated according to the distribution of materials on tracers55,57,59. Here, we further assume an incompressible viscoplastic rheology that is governed by the conservation of mass, moment, and energy (Eqs. (S1–S3) in the supplements). More details of the modeling method can be found in the supplements.

We presented a total number of 108 sensitivity tests to explore the SZI development. Following the initial conditions in previous studies11,42, we designed two types of numerical models for subduction evolution (with vs. without plate push, i.e., trenchward contraction Fig. 3), corresponding to the two scenarios in Fig. 1. In these models, the initial oceanic plate separates from a thin adjacent plate along a decoupling zone, e.g., a transform fault8,9,14,62,63. The key variables for the first type of model are the oceanic plate’s thermal age (80 or 100 Myr), the compositional density difference between the oceanic lithospheric mantle and the underlying asthenospheric mantle (by varying the asthenosphere’s compositional density, \({{{{{{\rm{\rho }}}}}}}_{{{{{{\rm{AM}}}}}}}\)), the minimum viscosity (\({{{{{{\rm{\eta }}}}}}}_{\min }\)) of the region following the initial oceanic plate, e.g., a mid-oceanic ridge or newly formed oceanic lithosphere (Fig. 6), and the slab’s maximum yielding stress (\({{{{{{\rm{\sigma }}}}}}}_{{{{{{\rm{slab}}}}}}}\)) (Supplementary Fig. S3). In the second type, the trenchward contraction is introduced by adding horizontal velocity at the end of the initial oceanic plate24,64 (Fig. 3b), which is active either throughout the modeling time or only during the first 2 Myr (Fig. 9); besides, \({{{{{{\rm{\eta }}}}}}}_{\min }\) and \({{{{{{\rm{\sigma }}}}}}}_{{{{{{\rm{slab}}}}}}}\) are also varied to test their effects on subduction continuity and average subduction rate (Fig. 9 and Supplementary Figs. S2 and S4).

In the first-type numerical models (Fig. 1a, b), where the slab lithospheric mantle shares the same compositional density with the underlying asthenospheric mantle (\({{{{{{\rm{\rho }}}}}}}_{{{{{{\rm{AM}}}}}}}={{{{{{\rm{\rho }}}}}}}_{{{{{{\rm{LM}}}}}}}\)), even though the ocean is old and cold (Figs. 3 and 6 and Supplementary Fig. S3), SZI cannot happen efficiently14,62. The gravity instability generally takes >6 Myr to grow enough to allow SZI (Fig. 6). In contrast, in runs where \({{{{{{\rm{\rho }}}}}}}_{{{{{{\rm{LM}}}}}}} > {{{{{{\rm{\rho }}}}}}}_{{{{{{\rm{AM}}}}}}}\), a less likely case based on mantle xenolith studies65, SZI occurs within 3 Myr (Fig. 6). Once SZI happens, the average subduction rate of the initial oceanic plate (\({{{{{{\rm{\upsilon }}}}}}}_{{{{{{\rm{ave}}}}}}}\)) is extremely high (i.e., >25 cm/yr), a universal feature of similar numerical models8. Only in 4 runs (4/48), \({{{{{{\rm{\upsilon }}}}}}}_{{{{{{\rm{ave}}}}}}}\) is <10 cm/yr (Fig. 6), i.e., the maximum rate of Pacific subduction after ~52 Ma49. Based on more sensitivity tests, by increasing ηmin, υave can be slowed down (Fig. 6). However, when ηmin is > 1022 Pa∙s, the subducted slab portion tends to break up with the following within 2–3 Myr after SZI24 (Fig. 6); subduction re-initiation cannot happen in the following >20 Myr because, after broad slab tearing, asthenospheric mantle upwells and melts suddenly, then solidifies to be a cold lithosphere, while the tip of the relict oceanic plate is thinned out and more buoyant than the new lithosphere (Fig. 10 and Supplementary Movie 2). Typically, the rifted upper plate, due to the SZI, is >300-km wide in these runs and keeps its width throughout the models11.

Fig. 10: A typical first-type model where slab tearing happens (Run 32, Fig. 6).
figure 10

a–h Model evolution regarding material, temperature, and velocity fields. A complete evolution of this model is in Supplementary Movie 2.

Full size image

By contrast, due to trenchward contraction, SZI happens within ~1.5 Myr in the second-type models (Figs. 1c, d). The subduction process is mostly continuous (22/32) except when \({{{{{{\rm{\eta }}}}}}}_{\min }\) is >1022.5 Pa∙s (Fig. 9), and in nearly half (15/32) of the runs, \({{{{{{\rm{\upsilon }}}}}}}_{{{{{{\rm{ave}}}}}}}\) can be <10 cm/yr (Fig. 9). Even if slab breakup occurs, subduction re-initiation can happen within 3 Myr in runs where trenchward contraction remains active24 (Supplementary Supplementary Movie 3). In this model type, the width of the upper-plate rift varies between 0 and ~300 km. Especially when the \({{{{{{\rm{\upsilon }}}}}}}_{{{{{{\rm{ave}}}}}}}\) is extremely high (>50 cm/yr), the rifting basin closes within 2 Myr after SZI due to the trench’s forward motion (Supplementary Fig. S5 and Supplementary Movie 4), in which case the melting and dehydration are limited14, and the forearc basin is erased23.

In addition, \({{{{{{\rm{\sigma }}}}}}}_{{{{{{\rm{slab}}}}}}}\) does not identifiably affect either model’s evolution (Supplementary Figs. S3 and S4). Regarding the second-type model, the duration and value of the horizontal velocity do not change the general pattern of model evolution, although they can mildly affect \({{{{{{\rm{\upsilon }}}}}}}_{{{{{{\rm{ave}}}}}}}\) (Supplementary Fig. S2). This is because, after ~225–300 km subduction, the slab’s down-going motion is dominated by its negative buoyancy53.