Research ArticleGEOPHYSICS

Lateral propagation–induced subduction initiation at passive continental margins controlled by preexisting lithospheric weakness

See allHide authors and affiliations

Science Advances  04 Mar 2020:
Vol. 6, no. 10, eaaz1048
DOI: 10.1126/sciadv.aaz1048


Understanding the conditions for forming new subduction zones at passive continental margins is important for understanding plate tectonics and the Wilson cycle. Previous models of subduction initiation (SI) at passive margins generally ignore effects due to the lateral transition from oceanic to continental lithosphere. Here, we use three-dimensional numerical models to study the possibility of propagating convergent plate margins from preexisting intraoceanic subduction zones along passive margins [subduction propagation (SP)]. Three possible regimes are achieved: (i) subducting slab tearing along a STEP fault, (ii) lateral propagation–induced SI at passive margin, and (iii) aborted SI with slab break-off. Passive margin SP requires a significant preexisting lithospheric weakness and a strong slab pull from neighboring subduction zones. The Atlantic passive margin to the north of Lesser Antilles could experience SP if it has a notable lithospheric weakness. In contrast, the Scotia subduction zone in the Southern Atlantic will most likely not propagate laterally.


It is well established that subduction is the driving mechanism for modern plate tectonics (15); however, the dynamics of subduction initiation (SI) is still poorly understood (68). New subduction zones may initiate in both intraoceanic and ocean-continent transitional (OCT) settings (9). The feasibility of SI at OCT regions (i.e. passive continental margins) is crucial for understanding the complete Wilson cycle (7). On the basis of theoretical and numerical analyses and natural observations (1012), previous studies have shown that it is difficult to initiate subduction in intact oceanic lithosphere and at passive continental margins. In addition, Nikolaeva et al. (13) and Marques et al. (14) used two-dimensional (2D) and 3D thermomechanical models, respectively, to investigate the possibility and controlling factors of SI at passive margins, which indicate that a thinned, hot continental lithosphere is required. However, the typical Atlantic-type passive margin is stable and, thus, difficult for SI. Recently, Baes and Sobolev (15) proposed that the SI at passive margins can be triggered by active mantle flow, originating from neighboring subduction zones and downgoing slab remnants from former subduction zones. These results indicate that unusual conditions are generally required for SI at passive margins due to the strength of typical OCT lithosphere (13, 16). Thus, understanding the transition from passive continental margins to subduction zones (active continental margins) remains as a major geoscientific challenge (7).

Intraoceanic SI is possible, either spontaneous or induced, only if sufficiently long lithospheric weak zones exist, e.g., transform faults or oceanic fracture zones (9, 12, 1719). Three intraoceanic subduction zones exist in the Atlantic: Scotia, Lesser Antilles, and Gibraltar (20, 21). These adjoin passive continental margins and these juxtapositions could provide favorable locations where existing plate convergence could propagate laterally to form new subduction zones (also called subduction invasion) (7, 22). However, whether SI at a typical passive margin can be induced by drag from a neighboring subduction zone remains uncertain. Possible cases of lateral propagation–induced SI [subduction propagation (SP)] in nature are summarized in Fig. 1. The Banda region (Indonesia subduction zones) illustrates a general way of SP (7, 23) (Fig. 1). The New Hebrides and Macquarie subduction zones may also be experiencing SP along weak zones (9, 24) (Fig. 1). Plate tectonic reconstruction suggests that SP accompanies Nazca plate subduction along the west coast of South America (25). Previous studies proposed that SP may also occur in the Atlantic passive margin (16, 2629); however, there is no evidence that SP is underway adjacent to the Scotia and Lesser Antilles convergent margins. Thus, the possibility of SP at the Atlantic passive margins remains unclear.

Fig. 1 Major convergent margins with possible lateral propagation of SI.

The preexisting faults can be transform faults or STEP. STEP, Subduction-Transform Edge Propagator. Topography data are from ETOPO1 Global Relief Model ( The topographic map is produced by the Generic Mapping Tools (64).

In this study, 3D thermomechanical models are conducted to investigate the possibility of SP at passive continental margins. The critical conditions include the strength of the OCT lithosphere, the age of the adjacent subducting slab, and the strain weakening of the lithosphere. Variable geodynamic regimes are manifested from the systematic numerical models, which are further compared to natural observations to identify regions of Atlantic passive margins where SP is likely.


Initial model setup

To quantify the possibility of SP at passive continental margins, a series of 3D numerical experiments are conducted with the initial model setup shown in Fig. 2. The initial model domain is 2000 km by 2000 km by 600 km, with spatial resolutions of 5.88, 5.88, and 2.46 km, respectively (x, y and z directions). Passive continental margins (OCT regions) are located on both sides of the convergent margin (Fig. 2C), with a preexisting intraoceanic subduction zone in between (Fig. 2B). An initial weak zone is prescribed for the intraoceanic convergent margin to decouple subducting and overriding oceanic plates. The passive continental margin is either stable (Fig. 2C) or with an initial 50-km-wide weak zone (Fig. 2D). Detailed model parameters and boundary conditions are shown in Materials and Methods and Supplementary Materials.

Fig. 2 Initial model setup illustrating geometry, lithologies, and thermal configurations.

(A) 3D model domain with colors indicating different lithologies as shown in the color bar at the bottom left. The top layer of sticky air is cut off for clarity. (B) XY slice at Z = 1000 km through the intraoceanic subduction zone. (C and D) XY slice at Z = 1750 km through the passive continental margins, with (D) or without (C) a preexisting lithospheric weak zone. The white lines in (B) to (D) denote isotherms, starting from 100° to 1300°C with 400°C intervals.

Model result with no weak zone at a passive continental margin

This model is conducted for passive continental margins without any preexisting OCT weak zones. Continued intraoceanic subduction is predicted with tearing on the edges of the subducting slab, i.e., the formation of STEP (Subduction-Transform Edge Propagator) fault (see Fig. 3) (30).

Fig. 3 Numerical model evolution with stable continental margin (no preexisting weak zone).

(A to C) Model snapshots at 2.55 Ma, 4.21 Ma, and 5.58 Ma, respectively. (Left column) Subducting slab morphology by temperature field. (Middle column) Relative surface elevation. (Right column) Second invariant of strain rate tensor of the horizontal depth slice at Y = 80 km. Tearing along the lateral edges of the subducting slab to form a STEP fault is predicted. See detailed model evolution in movie S1.

At the beginning, the subducting slab sinks deeper into the mantle under its own negative buoyancy (Fig. 3A). Subduction is accompanied by trench retreat and slab tearing adjacent to the neighboring continental lithosphere (Fig. 3, B and C). Models with variable ages and properties of the subducting lithosphere always predict similar STEP tearing, unless adjacent passive continental margins have weak zones. The strain rate map (the third column of Fig. 3) shows the STEP fault nucleated along the slab edge as a result of slab rollback. Formation of STEP faults along the edges of a subducting slab is a common phenomenon, with the Scotia convergent margin as an example (Fig. 1).

Model result with preexisting weak zones at a passive continental margin

Lithospheric weak zones such as transform faults, fracture zones, or rifts may exist along or near passive continental margins (3133), such as the continental transform fault on the southern margin of west Africa, the San Andreas fault system in California, and the Alpine Fault in New Zealand. To study the possible effects of these weak zones on SP, a numerical model with weakened passive continental margin was conducted (Fig. 4).

Fig. 4 Numerical model evolution with a weakened passive continental margin (preexisting weak zones in the OCT regions).

(A to D) Model snapshots at 2.81 Ma, 4.42 Ma, 5.79 Ma, and 6.96 Ma, respectively. For explanation of panel meanings, see Fig. 3 caption. Lateral propagation–induced SI at passive continental margin is predicted. See detailed model evolution in movie S2.

At the beginning, subduction in the central segment continues, with trench retreat (Fig. 4A). The trench retreat velocity is generally fast in this stage. With increasing slab pull, the subduction zone propagates laterally and leads to SI in neighboring passive continental margins (Fig. 4, B and C). After about 4 million years (Ma), subduction commences in the adjacent OCT regions (Fig. 4B). Because of contrasting convergence rates between the intraoceanic (fast) and OCT (slow) regions, a curved trench results (Fig. 4, B and C). As the model evolves, the subduction zone propagates further along the weak zone. After about 7 Ma, the whole plate begins to subduct (Fig. 4D). The rate of lateral propagation–induced SI at passive margins is about 17 cm/year, which is consistent with previous studies on intraoceanic SI propagation (22). The curved trench will gradually straighten as the entire slab subducts (cf. relative elevation in Fig. 4, C and D). According to the strain rate and velocity fields, subduction and trench retreat rates slow significantly during SI at passive margins because slab pull is partially applied to bend the adjacent region to initiate subduction.

Geodynamic regimes

SP is never predicted in existing numerical models for passive continental margins lacking a lithospheric weak zone (Fig. 2C). For this reason, we conducted a numerical experiment with such a weak zone (Fig. 2D). In addition, the effects of three main controlling factors are investigated, i.e., the age of subducting oceanic lithosphere, the friction coefficient of prescribed OCT weak zones at passive continental margins, and the strain weakening of lithospheric mantle. The numerical results (with OCT weak zones) are summarized as three contrasting modes (Fig. 5): subducting slab tearing along a STEP fault (STEP formation mode, green), lateral propagation–induced SI at passive continental margin (SI propagation mode, blue), and aborted subduction with slab break-off (slab break-off mode, red).

Fig. 5 Schematic regimes of 3D subduction mode selection.

(A to C) Regime diagrams correspond to different strain weakening effects. The effects of three controlling factors are investigated: the age of subducting oceanic lithosphere, the friction coefficient of passive continental margin weak zones, and the strain weakening of lithospheric mantle. Numerical results are summarized into three contrasting modes: subducting slab tearing along STEP faults [(D) green], lateral propagation–induced SI at passive continental margin [(E) blue], and aborted subduction with slab break-off [(F) red]. The typical model evolutions of (D) and (E) are shown in Figs. 3 and 4, respectively. The detailed model evolution of slab break-off (F) is shown in movie S3.

If strain weakening is not considered (Fig. 5A), then all three modes are predicted, depending on the age of oceanic lithosphere and the friction coefficient of the passive margin weak zone. The STEP formation mode is favored for old oceanic plates and stronger weak zones. Older oceanic lithosphere causes large slab pull, which drives subduction. However, stronger weak zones prevent induced SI at passive margins. Thus, the subducting slab tends to tear and decouple from a neighboring plate by forming STEP faults (Fig. 5D). In contrast, weaker zones at passive margins allow the large slab pull from older subducting slab to drag down the neighboring plate and propagate to form a new subduction zone along the OCT region (Fig. 5E). In cases with young subducting slab, slab pull force is low, which thus prevents either STEP fault formation (Fig. 5D) or lateral propagation (Fig. 5E). Instead, subduction will be aborted by break-off of the subducted young slab (Fig. 5F).

Strain weakening of lithospheric mantle is also important for mode selection (Fig. 5). With moderate strain weakening (Fig. 5B), 60 Ma old oceanic lithosphere will subduct continuously, leading to either STEP formation or SI propagation modes. However, the aborted subduction mode results from a nonzero friction coefficient for the OCT weak zone in the absence of strain weakening (cf. Fig. 5, A and B). With high strain weakening (Fig. 5C), the aborted subduction mode is never predicted. Younger oceanic slab and stronger weak zones lead to STEP fault formation mode. Alternatively, the lateral propagation–induced SI at passive continental margins results with older oceanic slab and lower strength of weak zones. In summary, strain weakening of lithospheric mantle contributes to strain localization and promotes subduction with STEP fault formation.

In the previous models, the two preexisting, symmetrically located OCT weak zones are always 50 km wide. To further study its importance, additional models with narrower OCT weak zones (i.e., 10 to 40 km) are conducted (fig. S1 and table S2). The results indicate that a narrow weak zone of 10 km prevents the SI propagation, differing from models with a 50-km-wide weak zone (fig. S1). However, the width of the OCT weak zone matters less when it is larger than 20 km, the model results of which are consistent with the reference models (fig. S1). Four additional asymmetrical models are further constructed with a prescribed single weak zone only in one of the neighboring passive margins (fig. S2 and table S2). The results indicate that SI propagation occurs in the OCT region with preexisting weak zone and old neighboring subducting slab. In contrast, the young subducting slab leads to aborted subduction with slab break-off. These results are, thus, consistent with the reference models in Fig. 5. On the other side without a weak zone, the passive margin is always stable without subduction collapse, leading to formation of a STEP fault. The prescribed initial length of the subducting slab is another factor that may influence what happens along the adjacent OCT region and is tested with additional models (fig. S3). The results indicate that the longer subducted slab leads to faster SP; however, it plays a minor role in determining whether or not subduction propagates into the OCT region (Fig. 5 and table S2).


This study explores the conditions whereby subduction can laterally propagate from an existing convergent margin into an adjacent passive continental margin. Numerical models reveal that a lithospheric weak zone is required for this to occur. In addition, an older and longer subducting oceanic slab is also important because this provides a greater slab pull force. The OCT region is a potential site for weak zones (34, 35), which could be caused by exhumed mantle serpentinite (36), inherited weakness (37), rifts, grain size damage, thermochemical plumes, and continental transforms (31, 32). These factors may result in a sufficiently weak OCT region to provide favorable conditions for SI to propagate if an active subduction zone is next to it. In a similar numerical study (38), SP and trench retreat were investigated with variable viscoplastic properties of the mantle, slab, and continental margin and were further applied to the tectonic evolution of the Western Mediterranean region. These models also indicated that strong slab pull in tandem with a weak continental margin resulted in SP, which is consistent with our models. Their model results also showed that the subducted slab teared with the strong continental margins. These may correspond to the STEP formation mode in the current study. However, the break-off mode is not popularly predicted in their study due to the old (~100 Ma) subducting oceanic plate.

SI propagation along a preexisting weak zone should be a common phenomenon of natural subduction systems (7, 22), e.g., the New Hebrides and Macquarie-Puysegur, as shown in Fig. 1. In the eastern part of the New Hebrides subduction system, the Matthew and Hunter subduction zones initiated around 2 Ma ago (6, 24, 39) and is generally regarded as an example of a transform fault evolving into a subduction zone driven by the neighboring sinking slab. Similarly, the existing active Puysegur subduction zone at the northern end of the Macquarie complex marks the boundary between the Australian and Pacific plates (40). Puysegur SI began within the last 10 to 12 Ma (41) and has propagated from north to south along the weak zone (evolving from spreading ridge to strike-slip fault and then to subduction zone) (9, 33, 42, 43).

Plate tectonic reconstructions show that the Banda subduction zone developed by the eastward propagation of the Java trench into an oceanic embayment by tearing along a former OCT boundary started at about 15 Ma ago (23, 44). An older subducted slab provided the driving force to drag down the neighboring oceanic lithosphere. Although there is no clear evidence for a preexisting weak zone in the OCT region, we propose, according to our numerical models, that the OCT region was weakened previously to accommodate SP into this region.

Our numerical model results may also have significant implications for SP in some other Pacific subduction systems, such as Cascades and Izu-Bonin-Mariana (Fig. 1). For the Cascades subduction system, the neighboring transform fault boundary favors SI propagation, although the subducting slab may be decoupled from neighboring oceanic plates by spreading ridges. Recent studies demonstrate that the Cascades subduction zone is propagating to the north along the Queen Charlotte margin (45, 46). The possibility of lateral propagation of the Cascades subduction zone requires further evaluation. The Yap subduction zone in the southern part of Izu-Bonin-Mariana initiated at about 19 to 15 Ma ago along a spreading ridge, which was a preexisting weak zone (12). Subduction failed later on the Yap convergent margin, which may be due to the blocking by attempted subduction of the thick and buoyant Caroline ridge (47). Note that some of the above natural cases are intraoceanic SI propagation, which are, thus, not exactly the same as what we have modeled.

Previous studies also propose that SP may occur in Atlantic passive margins from the existing Scotia and Lesser Antilles subduction zones (20, 48). Mart et al. (49) argued that two continent strips near the Lesser Antilles and Scotia may be weakened by the Pacific plate subduction-related melt and water infiltration. Duarte et al. (48) predicted that SP is likely near the Lesser-Antilles, Scotia, and southwest Iberia in the next 20 Ma. However, our results indicate that SP into adjacent passive margins is not as easy as thought previously. Subduction of the Atlantic plate beneath the Scotia plate began about 40 Ma ago (50, 51). Continuous trench retreat has formed STEP faults on both sides of the subducting slab (Fig. 1) (51, 52). According to our numerical models, the Scotia subduction zone is unlikely to propagate laterally to initiate a new subduction zone along the Argentine passive margin for two reasons: (i) the subducting slab is too young (<60 Ma), which does not favor SI propagation, and (ii) the formation of STEP faults prevents the lateral propagation of the existing subduction zone. Another active Atlantic subduction zone is the Lesser Antilles, for which SI occurred around 90 to 85 Ma ago (53). A continental transform margin now defines the northern limit of the Lesser Antilles convergent margin (32), which is a critical condition for SI propagation. If there is a weak zone on the eastern margin of the Bahama Platform, then the Lesser Antilles subduction zone may propagate to the north along the Atlantic passive margin. Teleseismic data show that thrust earthquakes are distributed along the Dominican region (northwest of the Lesser Antilles) to Cuba (fig. S4). These may indicate that subduction will someday propagate north of the Lesser Antilles subduction zone, first northwestward along the Dominican region, and then changing to northward into the Atlantic passive margins.

In summary, SI propagation along passive continental margins requires sufficiently weak OCT lithosphere and large slab pull force. Otherwise, STEP faulting or slab break-off will result (Fig. 5). After evaluating the main possible SI propagating cases in nature, we propose that SI propagating along a weak zone is the main way that a passive margin can be converted into a subduction zone. The Atlantic passive margin to the north of Lesser Antilles may or may not be vulnerable to SP, depending on the strength of the margin; this possibility needs further investigation. However, the Scotia subduction zone in the South Atlantic is unlikely to propagate due to the young subducting plate and the presence of a STEP fault along the edges of the subducting plate.


Model geometry and boundary conditions

The initial 3D model domain is 2000 km by 2000 km by 600 km (Fig. 2A) distributed over 341 × 341 × 245 nodes, with the spatial resolution of 5.88 km in the x direction, 5.88 km in the z direction, and 2.46 km in the y direction. Each grid cell has 2 × 2 × 2 randomly distributed markers used to trace rock composition, various material properties, and temperatures.

The velocity boundary conditions are free slip, except the permeable bottom boundary (5456). An external bottom boundary implies that the free-slip condition is satisfied at ~250 km below the base of the model domain. The external boundary condition allows global conservation of mass in the computational domain and is implemented by using the following limitations for velocities at the lower boundary of our models: ∂vx/∂y = 0, ∂vz/∂y = 0, ∂vy/∂y = (vyexternalvy)/∆yexternal, where ∆yexternal is the vertical distance from the lower boundary to the external boundary, and vyexternal = 0.

Compositional and thermal configurations

In the initial model setup (Fig. 2), the overriding plate is composed of an oceanic plate in the middle part and two adjacent continental lithospheres, which is similar to the Scotia subduction zone and Lesser Antilles subduction zone (Fig. 1). The subducted slab is an oceanic plate with variable ages. An initial weak zone is imposed between the subducting and overriding plates, which represents the preexisting subduction interface. The weak zone has a “wet olivine” rheology (57) in contrast to the “dry olivine” rheology elsewhere in both the lithospheric and asthenospheric mantle. For both the subducting and overriding oceanic plates, the initial material field consists of 3-km-thick upper oceanic crust, 5-km lower oceanic crust, and the oceanic lithosphere, the thickness of which depends on the half-space cooling model (58). The overriding continental plate incorporates 20-km upper continental crust, 15-km lower continental crust, and the underlying 63-km lithospheric mantle.

The top surface of the lithosphere is implemented by a layer of “sticky air” (5961), which is used to simulate an internal free surface. It is an initially 15-km-thick layer above the upper crust, which dynamically changes during experiments. The composition is either “air” (1 kg/m3, above y = −12 km water level) or “water” (1000 kg/m3, below y = −12 km water level), both of which have a low viscosity (1019 Pa s).

The thermal structure of the continental lithosphere varies linearly from 0° to 1300°C from the surface to the lithosphere bottom. The thermal structure of the oceanic lithosphere is applied using the half-space cooling model with given age (58). The initial temperature gradient in the asthenospheric mantle is about 0.5°C/km. The thermal boundary conditions have a fixed value (0°C) for the upper boundary and zero horizontal heat flux across the vertical boundaries. For the lower thermal boundary, an infinity-like external constant temperature condition is imposed, which allows both temperatures and vertical heat fluxes to vary along the permeable lower boundary of the box. This implies that the constant temperature condition is satisfied at ~250 km below the bottom of the model. This condition is implemented by using the equation Ty=(TexternalTy)/yexternal at the lower boundary, where Texternal is the temperature at the external boundary and ∆yexternal is the vertical distance from the lower boundary to the external boundary.

In this study, we use the I3ELVIS code (62) to solve the momentum, continuity, and energy equations. The code is based on the finite differences method and marker-in-cell techniques (62, 63).

Governing equations

The 3D Stokes equation for creeping flow is given by the followingσxxx+σxyy+σxzz=Px(1)σyxx+σyyy+σyzz=Pygρ(C,P,T)(2)σzxx+σzyy+∂σzzz=Pz(3)

Mass conservation is described by the incompressible continuity equationvxx+yy+vzz=0(4)

Heat conservation equations take the following formρCp(DTDt)=qxxqyyqzz+Hr+Ha+Hs(5)qx=kTx,qy=kTy,qz=kTz(6)Ha=TαDPDt(7)Hs=σxxε̇xx+σyyε̇yy+σzzε̇zz+2σxyε̇xy+2σxzε̇xz+2σyzε̇yz(8)where x, y, and z are the coordinates along the three directions, which are demonstrated in the initial model (Fig. 2); vx, vy, and vz are the velocity vectors; t is the time; σij are the viscous deviatoric stress tensor; ε̇ij are the components of the strain rate tensor; g is the gravitational acceleration, and the density ρ depends on composition (C), temperature (T), and pressure (P); k is the thermal conductivity; Cp is the isobaric heat capacity; and Hr, Ha, and Hs represent radioactive heat production, the adiabatic heating, and shear heating, respectively. For a given rock type, the density varies with pressure and temperature according to the relationρP,T=ρ0[1α(TT0)][1+β(PP0)](9)where ρ0 is the reference density at standard state P0 = 0.1 MPa and T0 = 298 K, and α = 3 × 10−5K−1 and β = 1 × 10−5MPa−1 are the thermal expansion and compressibility of all the rocks, respectively.

Rheological model

The viscoplastic flow is adopted in this study. The non-Newtonian viscous creep depends on the temperature, pressure, and strain rate (57). The ductile creep viscosity ηductile is defined by viscosity of dislocation (ηdisl) and diffusion (ηdiff) creep as1ηductile=1ηdisl+1ηdiff(10a)ηdisl=0.5(ε̇II)1nn(AD)1n exp(E+PVnRT)(10b)ηdiff=0.5(σcr)1nAD exp (E+PVRT)(10c)where P is the pressure, T is the temperature (in kelvin), and ε̇II is the second strain rate, which is defined as ε̇II=(12ε̇ijε̇ij)1/2. AD is the pre-exponential factor, E is the activation energy, V is the activation volume, and n is the creep exponent. These parameters are determined by the flow law experiments (table S1). σcr = 30,000 Pa is the critical stress for transition from diffusion to dislocation creep. R is the gas constant.

Besides viscous creep, the extended Drucker-Prager yield criterion is adopted to describe the plastic deformation mechanismσyield=C0+Pφ(11)ηplastic=σyield2ε̇II(12)where, σyield denotes the yield stress, C0 is the cohesion, P is the dynamic pressure, and φ is the friction coefficient of dry rocks. The fracture-related strain weakening of the mantle is implemented, with varied strain weakening parameters shown in table S2.φ={φa+(φbφa)×γγ0,if γ  γ0φb,if γ > γ0(13)andγ=12(ε·ij(plastic))2dt(14)where γ is the integrated plastic strain, γ0 is the upper limit strain weakening, and φa and φb are the initial and final friction coefficient values, respectively.

Last, the ductile rheology is combined with the brittle/plastic rheology to achieve an effective viscoplastic rheology.ηeff=min(ηductile,ηplastic)(15)

The lower and upper viscosity cutoff values are 1018 and 1024 for all kinds of rocks. The material properties and used parameters in this study are shown in tables S1 and S2.


Supplementary material for this article is available at

Fig. S1. Additional model results with variable widths (10 to 40 km) of the OCT weak zones, in comparison to the reference models with preexisting OCT weak zones of 50 km (i.e., top).

Fig. S2. Additional models with a prescribed weak zone only existing in one of the neighboring passive margins.

Fig. S3. Additional models with variable lengths of the preexisting subducted oceanic slab.

Fig. S4. The tectonic map of the Caribbean region.

Table S1. Physical properties of materials used in the numerical experiments.

Table S2. Parameters of the numerical models and the results.

Movie S1. Animation of the STEP formation mode.

Movie S2. Animation of the SI propagation mode.

Movie S3. Animation of the aborted SI with slab break-off mode.

References (6567)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Thorough and constructive reviews by M.-A. Gutscher and an anonymous referee improved the paper. Funding: This work was supported by the NSFC Tethys project (91855208), the Strategic Priority Research Program (B) of the Chinese Academy of Sciences (XDB18000000), the NSFC Excellent Young Scholar Project (41622404), and the support from the China Scholarship Council and University of Chinese Academy of Sciences. T.V.G. acknowledges support from EU projects CREEP and SUBITOP. This is UTD Geosciences contribution number 1352. Numerical simulations were run with the clusters of the National Supercomputer Center in Guangzhou (Tianhe-II). Author contributions: Z.-H.L. conceived the study. X.Z. and Z.-H.L. designed the numerical experiments. X.Z. performed the experiments. X.Z and Z.-H.L. analyzed the data. X.Z., Z.-H.L, T.V.G., and R.J.S. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article