Strike slip motion and the triggering of subduction initiation

Plate tectonic reconstructions of three of the best-defined Cenozoic subduction initiation (SI) events in the western Pacific, Izu-Bonin-Mariana, Vanuatu, and Puysegur subduction zones, show substantial components of strike-slip motion before and during the subduction initiation. Using computational models, we show that strike-slip motion has a large influence on the effective strength of incipient margins and the ease of subduction initiation. The parameter space associated with visco-elasto-plastic rheologies, plate weakening, and plate forces and kinematics is explored and we show that subduction initiates more easily with a higher force, a faster weakening, or greater strike-slip motion. With the analytical solution, we demonstrate that the effect of strike-slip motion can be equivalently represented by a modified weakening rate. Along transpressive margins, we show that a block of oceanic crust can become trapped between a new thrust fault and the antecedent strike-slip fault and is consistent with structural reconstructions and gravity models of the Puysegur margin. Together, models and observations suggest that subduction initiation can be triggered when margins become progressively weakened to the point that the resisting forces become smaller than the driving forces, and as the negative buoyancy builds up, the intraplate stress eventually turns from compressional into extensional. The analytical formulation of the initiation time, t_SI, marking the moment when intraplate stress flips sign, is validated with a computational models. The analytical solution shows that t_SI is dominated by convergence velocity, while the plate age, strike-slip velocity, and weakening rate all have a smaller but still important effect on the time scale of subduction initiation.

1 Introduction

Consensus on a unified description of subduction initiation has been slow to develop as initiation is a transient process whose record is generally obscured by subsequent subduction zone processes, notably burial, overprinting, uplift, and compression and over-thrusting. Nevertheless, there is a substantial geological record with nearly all ocean-ocean subduction zones having initiated since the end of the Mesozoic and about half of all ocean-continent ones having re-initiated since the mid-Mesozoic (Hu and Gurnis, 2020). Subduction initiation, moreover, occurs nearby existing subduction zones (Crameri et al., 2020), is a fundamental component of plate tectonics, and is putatively associated with key changes in the force balance of tectonic plates. For example, the Pacific Plate changed its direction of motion from NNW to NW at around 50 Ma (Whittaker et al., 2007; Torsvik et al., 2017), synchronously with the initiation of two major subduction zones in the western Pacific, the Izu-Bonin-Mariana (IBM) (Ishizuka et al., 2018; Reagan et al., 2019) and Tonga-Kermadec (Sutherland et al., 2020). However, why new subduction zones form remains poorly understood, as subduction initiation appears to be mechanically unfavorable with initial slab pull being insufficient to overcome resistance from the friction between plates and bending of the slab (McKenzie, 1977; Toth and Gurnis, 1998; Li and Gurnis, 2022). An external force, or low initial strength between plates, is required to start subduction, with scenarios for initiation described as either spontaneous (Stern and Bloomer, 1992; Nikolaeva et al., 2010) (with no external force but low strength) or induced (Toth and Gurnis, 1998; Gurnis et al., 2004) (with an external force).

With theoretical and computational approaches, the thermal age of plates, compositional variations, trench-normal convergence and in–plane stress (Gurnis et al., 2004; Nikolaeva et al., 2010; Leng and Gurnis, 2011) and fault strength, fault weakening and plate bending (McKenzie, 1977; Toth and Gurnis, 1998; Thielmann and Kaus, 2012; Qing et al., 2021) having been identified as key mechanical factors which respectively drive and limit subduction initiation. Independent of spontaneous and induced scenarios, driving forces must overcome frictional resistance between plates and bending of the high, effective viscosity plate. If the plate boundary does not weaken sufficiently fast, the oceanic plate will not slide into the mantle to allow the negative thermal buoyancy to grow quickly enough. Constitutive models with strength that weakens with deformation, due to grain size reduction, grain damage, or volatile ingestion, have been identified (Hirth and Kohlstedt, 1996; Thielmann and Kaus, 2012; Bercovici and Ricard, 2014) and can lead to rapid instability (Leng and Gurnis, 2015; Zhou et al., 2018). Models exhibiting spontaneous initiation start in a critical state that verges on instability (Leng and Gurnis, 2015; Zhou et al., 2018), an unsatisfactory condition to address causes for the onset of initiation. If a plate boundary is in a state close to instability, the question arises as to why the boundary initiated at that point in time and not earlier.

Strike-slip motion might be one reason plate boundaries are brought closer to a state favorable for subduction initiation. Plate tectonic reconstructions of the IBM, Vanuatu and Puysegur subduction zones, three of the best constrained subduction initiation events in the western Pacific during the Cenozoic, each show a substantial component of strike-slip motion before and during initiation (Figure 1). The Puysegur subduction zone extending south from the South Island of New Zealand has a well-documented component of strike-slip motion during subduction initiation which continues to the present (Lamarche and Lebrun, 2000; Sutherland et al., 2006). The initiation of IBM, among the best studied and often used as a point of comparison with other subduction zones, ophiolite–origin models, and mechanical models (Arculus et al., 2019). Although less discussed, IBM experienced a strong component of strike slip motion during the well-documented period of initiation (Gurnis, 2023). The Vanuatu subduction initiation also saw a strong component of strike-slip motion during the interval 15 to 12 Ma when the new subduction zone was forming through a polarity reversal with velocities of about 5–6 cm/yr (Figures 1B, D). In addition to the main Vanuatu subduction zone, a new segment of the plate boundary is initiating at its southern boundary, referred to as the Matthew and Hunter subduction zone, with a small 2 cm/yr convergence and a subtantially larger strike-slip motion along the Hunter Ridge since 2 Ma (Patriat et al., 2019). All of these subduction initiation events must have been shaped by strike-slip motion, but do large components of strike-slip motion influence the mechanics of initiation?

FIGURE 1

FIGURE 1. Three Cenozoic subduction initiation events that occurred in the western Pacific showing a substantial component of relative, strike-slip motions at the nascent trench just prior to initiation. (A). Present-day western Pacific with the structures against which subduction initiated shown in brown. Regions of detail shown with boxes. (B). Magnitude of strike slip motion (above) and convergence direction (below, with 0 being normal to the strike of the trench) before onset of initiation. Izu-Bonin-Mariana (IBM) shown at 55 Ma (red) and 50 Ma (black solid) with values for intermediate times shown with grey shading, Vanuatu at 15 Ma (blue dashed) and Puysegur (black, dot-dashed) using the plate model of Müller et al. (2019). (C). Conditions at IBM at 55 Ma, 5 Myr before the onset of SI. (D). Vanuatu at 15 Ma, 1–3 Myr before onset. (E). Puysegur at 15 Ma, 1–3 Myr before onset. Abbreviations for plates: PA, Pacific Plate; AU, Australia Plate; PS, Philippine Sea Plate.

2 Model formulation

We solve for the Stokes equations using traditional formulations used in geodynamics (Moresi et al., 2000; Ismail-Zadeh and Tackley, 2010) and show that a component of strike-slip motion can substantially reduce the strength of a nascent plate boundary while providing a triggering mechanism to nucleate a new subduction zone. The idea that strike-slip motion produces a more favorable condition for the far-field compression to induce subduction is examined by Zhong and Li (2023) with 3D models. Here, we address this problem with a sliced 3D geometry extended from the trench-perpendicular cross-section (x-z dimension), with an additional trench-parallel dimension (y-dimension) that accounts for the strike-slip motion (Supplementary Figure S1). In the trench-perpendicular (x-z) dimension, we apply either a convergent velocity or a convergent force on the right side of the subducting plate to induce subduction initiation. In the trench-parallel dimension, a strike-slip velocity is applied on the right wall of the subducting plate to drive the strike-slip velocity. For the two boundaries normal to the strike direction, we apply a periodic boundary condition, which allows the material to flow through the two boundaries freely. As the width of the trench parallel dimension is small ($\approx 6$ km) with only two layers of elements and periodic boundary conditions enforcing the strike-slip velocity, the models virtually solve for a plane strain problem with no variation along strike.

A younger overriding plate is to the left of an old oceanic plate with either a compressional force or a convergent velocity added on the right edge of the older plate. A pre-existing weak zone with yielding stress being reduced by half is located at the plate boundary. As plate convergence and strike-slip motion accumulates, the strain will eventually localize at the plate boundary.

The deformation of a visco-elasto-plastic material is computed with the finite element code Underworld (Mansour et al., 2019). The constitutive relationship between strain rate, $\dot{\epsilon }$, and deviatoric stress, τ, is defined through a visco-elastic Maxwell body (Moresi et al., 2003),

$\dot{\epsilon }=\frac{\dot{\tau }}{2\mu }+\frac{\tau }{2\eta }(1)$

Viscosity follows the non-Newtonian Arrhenius law,

$\eta ={\eta }_0{\left(\frac{{\dot{\epsilon }}_{II}}{{\dot{\epsilon }}_0}\right)}^{\frac{1}{n}-1}{e}^{\frac{E}{nR}\left(\frac{1}{T}-\frac{1}{T_0}\right)}(2)$

where μ is the shear modulus, η the viscosity, T the temperature, and ${\dot{\epsilon }}_{II}$ the second invariant of strain rate tensor. E, n and R are activation energy, non-linear exponent, and ideal gas constant. ${\dot{\epsilon }}_0$, η₀ and T₀ are reference strain rate, reference viscosity and reference temperature.

Plasticity, describing the strength of the rock, is an essential component of subduction initiation. We assume a yielding envelope given by the Drucker-Prager failure criteria bounded by a maximum stress.

$\tau =\min \left(C\cos \varphi +P\sin \varphi ,{\tau }_{max}\right)(3)$

where τ is the yielding envelop, C and ϕ are cohesion and friction angle and τ_max the maximum stress the rock can sustain. With the accumulation of plastic strain, the yielding envelope is reduced through rock damage or grain size reduction. The weakening is approximated as a two-stage process: The yielding envelope first reduces linearly with increasing plastic strain until saturation at which point the plastic parameters C and ϕ remain constant.

$C=C_0+\left(C_f-C_0\right)\min \left(1,\frac{{\epsilon }_P}{{\epsilon }_{P0}}\right)(4)$

$\varphi ={\varphi }_0+\left({\varphi }_f-{\varphi }_0\right)\min \left(1,\frac{{\epsilon }_P}{{\epsilon }_{P0}}\right)(5)$

The reference plastic strain, ɛ_P0, controls the weakening rate with rocks weakening faster with smaller ɛ_P0 and vice versa. With a lack of consensus on the underlying mechanisms responsible for weakening (Hirth and Kohlstedt, 1996; Thielmann and Kaus, 2012; Bercovici and Ricard, 2014), ɛ_P0 is varied as an unknown along with other key parameters including the external compression, F_xx, and strike-slip velocity, V_SS. Details of model parameters are in Supplementary Table S1.

In addition to the numerical models, we formulate an analytical solution modified from Li and Gurnis (2022) considering the extra influence of the strike-slip motions. The formulation of the analytical solution is based on the horizontal force balance of the subducting plate between the forces that drive the plate motion and forces that resists plate motion. The strike-slip motion affects the force balance equation through accelerating the weakening process, thereby reducing the frictional resistance at the plate boundary. By varying modeling parameters of interest, like the strike-slip velocity and weakening rate, we measure quantities like the plate force, F_xx, initiation time, t_SI, and total work, W_SI, from numerical models, that are compared with the analytical predictions.

3 Results

3.1 Kinematic boundary condition

We use a kinematic boundary condition with a constant convergent velocity on the right end of the subducting plate to initiate subduction (Supplementary Figure S1A; Figure 2A). The different conditions controlling the evolution of subduction are evaluated by tracking the vertically–integrated horizontal compressional stress in the plate, i.e., $F_{xx}={\int }_0^{d_{\mathit{lith}}}-{\sigma }_{xx}dz$, where d_lith is the thickness of lithosphere. Prior to initiation, plate motion is resisted by a large coupling stress at the plate boundary. Later, F_xx drops with plate convergence, ℓ, through plastic weakening and decoupling of the two plates and accumulation of negative buoyancy (Li and Gurnis, 2022). Eventually, F_xx becomes negative (extensional), indicating that subduction has become self-sustaining and driven by the negative buoyancy of the slab instead of external forces (Figures 2B, C), and we define the time when F_xx drops to 0 as the initiation time t_SI. Following Li and Gurnis (2022), we define the total work, W_SI, done by the boundary velocity to induce a subduction initiation until t_SI (i.e., $W_{SI}={\int }_0^{t_{SI}}{F}_{xx}{V}_{x}dt$), characterizing the total resistance that the driving force overcomes. The tracking of plate force, initiation time, and total work quantitatively reflect the difficulty to induce subduction initiation. A larger F_xx, t_SI and W_SI indicate the subduction initiation encounters greater resistance. We develop an analytical solution of the strike-slip subduction initiation model, validate it with the numerical results, and expand the parametric space via the analytical solution.

FIGURE 2

FIGURE 2. (A) An example of model with 1 cm/yr convergent velocity. Color shows the effective viscosity. Vectors show the in-plane component velocity and contours show strike slip velocity in cm/yr. (B) The integrated horizontal compressional force in the plate F_xx as a function of plate motion ℓ for models with different strike slip velocity V_SS under ɛ_P0 = 1. Solid line for numerical model results and dashed line for analytical results. (C) Same as (B) but with slower weakening rate ɛ_P0 = 2. (D) The work done to initiate a subduction W_SI of every numerical model (color-coded solid circle) in the parametric grid search varying the strike slip velocity V_SS and ɛ_P0. Contours show the prediction of W_SI from analytical model.

In 2D cross-sectional models (with no strike-slip motion), F_xx can be formulated analytically via a force balance analysis on the subducting plate, such that F_xx, together with slab pull as the driving force, is balanced by the resistance from plate bending, inter-plate friction, isostatic gradients at the boundary, and shear from mantle flow (Li and Gurnis, 2022). The addition of the third dimension (strike-slip) mainly changes the system by enhancing the weakening process. To analytically address the role of strike slip motion, we modify the formulation of inter-plate friction by enhancing the weakening as a result of the strike-slip velocity (see Supplementary Material S1). The contribution of strike-slip motion is analytically shown to be equivalent to reducing the effective weakening rate ${\epsilon }_{P0\text{\_}eff}={\epsilon }_{P0}\frac{V_x}{\sqrt{{V_x}^2+{V_{SS}}^2}}$ and thereby reduce the inter-plate friction. This friction is balanced by the driving force, F_xx, and the analytical model (dashed line in Figures 2B, C) predicts a decrease of F_xx with increasing V_SS, in agreement with the numerical models (solid line in Figures 2B, C). Increasing strike-slip velocity and weakening rate consistently lowers F_xx, as both reduce the resistance inhibiting subduction initiation. A grid search of V_SS and ɛ_P0 shows how these parameters influence W_SI (Figure 2D). In the limit of a small V_SS and a large ɛ_P0, convergence is insufficient to nucleate a weak zone and initiate a subduction (open circles, Figure 2D). Outside of that domain, subduction initiates and the resistance has a consistent trend: Increased V_SS and decreased ɛ_P0 tend to lower W_SI.

Finally, we evaluate t_SI from the analytical model varying V_x, V_SS, ɛ_P0 and subducting plate age (Figure 3). As discussed in Li and Gurnis (2022), the initiation time is dominantly controlled by the total plate convergence, so that convergence rate V_x has a first-order control on the predicted t_SI. Subducting plate age determines the plate thickness, which controls both the plate bending resistance and the slab pull. The plate age affects the force balance in two ways, the slab pull and plate bending because the plate thickness changes. As the dependence of slab pull on plate thickness is first order while plate rigidity (with plasticity) second order, the change of t_SI is dominated by slab pull when plate is young ($<10$ Ma), and for an old plate ($>20$ Ma) the t_SI is governed by plate bending. Therefore, below 10 Ma we see a decrease of t_SI as the plate gets older due to increased slab pull that drives plate motion, while above 20 Ma t_SI increases with plate age due to increased bending resistance. Both V_SS and ɛ_P0 influence t_SI via the term of inter-plate friction. Both increase of V_SS and decrease of ɛ_P0 lead to an acceleration of subduction initiation, thereby reducing the t_SI. However, compared to V_x, the impact of the plate age, V_SS and ɛ_P0 is minor.

FIGURE 3

FIGURE 3. Analytical prediction of initiation time t_SI under different convergent velocity Vx and subducting plate age with (A). V_SS = 0, ɛ_P0 = 1, (B). V_SS = 3 cm/yr, (C). V_SS = 6 cm/yr, ɛ_P0 = 1. (D–F) are the same as (A–C) except ɛ_P0 = 2. Red, magenta, blue and green box show the parameter range of IBM (50 Ma), IBM (55 Ma), Puysegur (PU) and Vanuatu (VT) subduction initiation from the given V_SS with a ±1 cm/yr range.

3.2 Force boundary conditions

Besides an applied convergent velocity, subduction initiation can be induced through applying a force on the edge of the plate (Supplementary Figure S1B). We use a geometry approximating a mid-ocean ridge and a compressional force F_xx at the ridge edge of the subducting plate (Supplementary Figure S1). The total compression driving the plate motion is composed of the internal ridge push from the thermal contrast and the external force F_xx. When F_xx is sufficiently large, plate bending and resistance by shearing between plates at the nucleating boundary are overcome with subduction initiation (Leng and Gurnis, 2011; Zhong and Li, 2019); the minimum force needed for initiation is determined through systematic variation of parameters.

All prior models of subduction initiation in the literature have no strike-slip motion, and so we start with a case with no strike-slip velocity with (ɛ_P0, V_SS, f_xx) = (2, 0 cm/yr, 8.57 × 10¹² N/m), as a reference (Supplementary Figure S2B) In this case, due to the insufficient compression to induce subduction initiation (F_xx = 8.57 × 10¹² N/m) under a relatively low weakening rate (ɛ_P0 = 2, larger ɛ_P0 meaning slower weakening), the subduction initiation fails to occur and the plate boundary remains stable (Supplementary Figure S2B). Typical ridge push forces range from 2 to 4 × 10¹² N/m (Bott, 1991; Toth and Gurnis, 1998) and so the F_xx used is not particularly small. Over time, the plate boundary remains stable with no subduction being induced. By taking this case as a reference, and then by either adding a V_ss of 2 cm/yr (Supplementary Figure S2C), or decreasing ɛ_P0 from 2 to 1 (such that the fault weakens faster), the plate boundary initiates into a subduction zone (Supplementary Figure S2D). This shows that the presence of strike-slip motion and acceleration of the weakening process independently facilitate subduction initiation.

The key parameters, including strike-slip velocity, external compression, and the weakening rate, are systematically varied while computing the initiation time (Supplementary Figure S3). The grid search shows that regardless of the choice of weakening rate and resolution, increased strike-slip velocities lower the minimum required compressional force (Supplementary Figure S3 thick lines) for subduction to initiate. A lower ɛ_P0 (faster weakening) always reduces the required external force, therefore facilitating subduction initiation. Unlike the velocity boundary model where the strike-slip velocity has a minor influence on plate stress and initiation time, the addition of strike-slip motion significantly reduces the stress required to initiate a subduction zone (solid curves in Supplementary Figures S3, S4), indicating the strike-slip motion provide a more favorable condition for subduction initiation.

3.3 Fault evolution

The computations yield a fault system that evolves from a single strike-slip fault into a strain partitioned system composed of a vertical strike-slip fault and a dipping oblique thrust fault with partitioning of strike-slip velocity between the two (Figure 4). Initially, the plate boundary fault is vertical, and will remain so with only subsequent strike-slip motion; however, subduction requires a dipping thrust fault to decouple the two plates. In those cases where the system evolves into a subduction zone, we have found that the thrust fault emerges through the combined reuse of the strike-slip fault and new fault formation. The upper part of a trapped block is bound by a vertical strike-slip fault and an oblique thrust fault as the upper lithosphere has a high viscosity and brittle failure (faulting) is the dominant mechanism that accommodates the compression. The lower part of the vertical strike-slip fault rotates with compression of the somewhat lower viscosity, ductile lower lithosphere. Through this process, the lower part of the vertical strike-slip fault is transferred and becomes the lower part of the oblique thrust fault. The reuse of a vertical fault weakened by strike-slip motion is a primary reason why strike-slip motion can lower the required force to induce subduction and therefore facilitate subduction initiation.

FIGURE 4

FIGURE 4. Fault evolution of a high resolution models with (ɛ_P0, V_SS, V_p) = (2, 6 cm/yr, 4 cm/yr). (A–D) Red and blue for oceanic and continental crust. Background color for temperature. (E–H) Second invariant of strain rate $({\dot{\epsilon }}_{II})$.

In the case shown above, we choose a high yield stress (300 MPa) which makes the crust of the overriding plate strong. However, with a continental upper plate with a weaker crustal layer, the resultant morphology is different (Supplementary Figure S5). In this case, the rheology yields a typical continental strength envelope (Kohlstedt et al., 1995), with a weak layer present at the base of the continental crust, decoupling the lithosphere from the crust. Unlike the strong plate case where a new dipping fault emerges in the upper part of the subducting plate in response to the far–field compression, in the case of weak continental crust, the overriding plate severely deforms to accommodate the compression. In the overriding plate, the lower lithosphere delaminates at the base of the buoyant continental crust. The compression forces the subducting plate and the lower lithosphere of the upper plate to bend simultaneously, and a vertical strike-slip fault emerges in the upper plate crustal layer to reconcile the strike-slip motion. Consequently, a strain partitioning system emerges, but the wedge trapped between the strike-slip fault and the dipping thrust fault is composed of continental crust.

4 Discussion

With the numerical models that extend a traditional 2D trench perpendicular cross-sectional domain with orthogonal strike–slip motion, such motion influences subduction initiation. This builds on the concept that mature strike-slip faults like the San Andreas Fault (SAF) are weak [as revealed through stress indicators (Zoback et al., 1987) and heat flow (Brune et al., 1969)] and that the damage around faults grows with increasing fault displacement (Faulkner et al., 2011; Savage and Brodsky, 2011). A common conceptual idea is that a strike-slip fault will have local irregularities along the strike such there will be zones of convergence that could be sites for subduction initiation. A well-known example where local transpression leads to downwelling, although not subduction initiation, is the Miocene evolution of the San Andreas Fault (SAF) and Transverse Ranges of southern California in North America. Here, it is thought that the clear convergence across the Big Bend segment of the SAF north of Los Angeles has led to crustal and lithospheric thickening and convective instability (Humphreys and Clayton, 1990). However, we present a scenario that is distinctly different from this one, advancing the hypothesis that even with no irregularities in the orientation of faults along strike, subduction initiation is enhanced by strike-slip motion.

Can the models be applied to specific subduction initiation events? An obvious application is along the boundary between the Pacific and Australian Plates, south of New Zealand, the Puysegur subduction zone (Figure 1E). Here, a major strike-slip fault immediately to the east of the trench accommodates some of the relative plate motion (Collot et al., 1995; Lamarche and Lebrun, 2000). The age of the Puysegur subduction zone has been estimated in two ways. First is through a reconstruction using a combination of the depth extent of the seismic slab and plate kinematics, which date the onset of subduction to between 15 and 12 Ma (Sutherland et al., 2006). Second, the age has been estimated through detailed seismic imaging of the stratigraphy on the overriding Pacific Plate which shows basin inversion (Shuck et al., 2022). In the northern section of Puysegur, the compression starts at 15 Ma and transitions to mild extension by 8 Ma. During this interval, most of the relative motion between the Australian and Pacific Plates was strike-slip with velocities of about 3 cm/yr (Figure 1B). The subducting plate at Puysegur Trench formed along the Macquarie Ridge Complex (MRC) system, where sea–floor spreading was active since the Eocene (about 40 Ma) but terminating at about 15 Ma (Lebrun et al., 2003). The subducting plate age ranges from 25 Ma (north) to 5 Ma (south). With the given range of V_x, V_SS, and plate age for Puysegur, we predict a range of t_SI from the analytical model (Figure 3, blue box). The northern Puysegur trench corresponds to the top right corner of the blue box (V_x = 1 cm/yr, plate age = 25 Ma), yielding a model t_SI ≈ 12 Myr which roughly agrees with the observed 8 Myr duration for the basin inversion (e.g., 16 to 8 Ma). For the southern Puysegur, the analytical model yields a larger t_SI ($>30$ Myr) given the small V_x (≈0.2 cm/yr) and young subducting plate (5 Myr old) at 15 Ma. This prediction may over-estimate t_SI as during the subduction initiation the convergent velocity V_x increases substantially from 0.2 cm/yr at 15 Ma to $\sim 2$ cm/yr today due to the change of plate motion direction (Choi et al., 2017), and the large predicted t_SI sheds light on why the southern Puysegur remains under compression today—in addition to the hypothesis of southward propagation of subduction initiation proposed in Shuck et al. (2022).

Puysegur has other features reflective of the initiation phase as it evolved from a strike-slip boundary to transpressional since 16 Ma, in particular a sliver of oceanic crust (Hightower et al., 2020) trapped at the Puysegur Trench, bounded by a vertical fault and an oblique thrust fault (Collot et al., 1995; Shuck et al., 2022). East of the Puysegur Trench, the morphology of the Puysegur Fault through the central part of the southern Puysegur Ridge is sharp and consistent with present-day activity (Shuck et al., 2021). The 2009 Mw 7.8 Dusky Sound earthquake below Fiordland, the strike-slip onshore extension of the nascent subduction zone, showed that nearly all of the relative motion, including the strike-slip motion, was accommodated on the thrust interface (Beavan et al., 2010). Consequently, the seismic imaging and 2009 event for Puysegur shows that it has a structure, morphology, and strain partitioning consistent with that shown in Figure 4, that is with substantial strike-slip motion on the thrust interface, and the entrapment of an oceanic sliver between the newly formed thrust interface and the strike-slip fault.

Subduction along IBM initiated at the boundary between the Pacific Plate and the smaller West Philippine Sea Plate along the Kyushu Palau Ridge (Figure 1C). The new subduction zone is inferred to have formed between 52 and 50 Ma as indicated by ages within the forearc (Reagan et al., 2019) and basement of the West Philippine Sea Plate (Ishizuka et al., 2018). IBM is often viewed as a type locality for studying subduction initiation because of the rock record on the IBM fore–arc and an association with changes in Pacific Plate motion during the Cenozoic (Stern and Bloomer, 1992; Arculus et al., 2019). The IBM subduction zone may be an example of spontaneous subduction initiation (Stern and Bloomer, 1992), with a short, rapid burst of localized extension (Reagan et al., 2019). The overriding Philippine Sea Plate had a relic arc with thick buoyant crust adjacent to the old, cold Pacific plate at the location of apparent subduction initiation and this may have played a role in its initiation (Leng and Gurnis, 2015). Computational models suggest that compression of the incipient plate margin needs to overcome a large resisting force associated with bending the incipient slab (Gurnis et al., 2004; Li and Gurnis, 2022), but there is debate in the literature on the significance of the compression. On the one hand, Maunder et al. (2020) argued that early compression cannot lead to near–field extension and boninite formation and that a large vertical force is needed for such extension, a force which only can exsit if there is pre-existing subduction. Moreover, a recent interpretation of Baron isotope systematics on IBM forearc recovered samples apparently indicate an initial low–angle phase of thrusting prior to the formation of the distinctive basalt to boninite sequence (Li et al., 2022).

Although not widely discussed, several lines of evidence suggest that there was a strong component of strike-slip motion during IBM inception. Paleomagnetic observations show that since 50 Ma, the Philippine Sea Plate experienced nearly 90° of clockwise rotation (Hall et al., 1995); since the orientation of the Kyushu-Palau Ridge (the relic arc that formed at the new boundary) is currently N-S, it would have been more E-W during subduction initiation. The orientation of convergence between the West Philippine Sea and Pacific Plates depends on the absolute motion of the Pacific Plate, but the Pacific moved by the overriding plate in a mostly strike-slip orientation with velocities of about 5 cm/yr for at least five million years before initiation (Figures 1B, C), with V_SS and V_x varying substantially along strike due to the curved plate boundary. Independently, detailed bathymetry just to the west of the Kyushu Palau Ridge reveal a strike–slip fault subparallel and normal faults perpendicular to the ridge in which the basement dates to the time of subduction inception, 49 Ma (Gurnis, 2023). This is consistent with transverse motion along the KPR at the time of initiation. Nearly identical forearc volcanic stratigraphy near the Bonin Islands and near Gaum Islands, 500 km apart when plate tectonic motions are accounted for (Leng and Gurnis, 2015), could have been further seperated by strike slip motion.

The strike-slip model of induced subduction initiation can be applied to IBM. With subducting plate age varying between 40 and 60 Myr old at 50 Ma (Hall et al., 2003; Lallemand and Arcay, 2021), the analytical model yields a range of initiation times (Figure 3). At 55 Ma, the predicted t_SI range from 5 to 10 Myr with V_SS =3–6 cm/yr. At 50 Ma, the reconstruction gives a larger range of V_SS (0–6 cm/yr), and the analytical model predicts a t_SI ≈ 3–5 Myr. Using the plate tectonic constraints at 55 and 50 Ma, we estimate the time for the IBM to experience conversion from compression to extension to be around 50 to 45 Ma, a 5 Myr duration slightly larger than the 2 Myr inferred from fore–arc opening determined from ⁴⁰Ar/³⁹Ar and U-Pb ages of the boninite-basalt sequence in the IBM fore–arc (52–50 Ma) (Reagan et al., 2019).

IBM subduction initiation may be associated with the shift in Pacific Plate motion from more northerly before 50 Ma to more westerly after 50 Ma, as evident from plate reconstructions and the bend in the Hawaiian–Emperor Seamount chain (e.g., Torsvik et al., 2017; Müller et al., 2019). Hu et al. (2022) showed that Pacific Plate motion could change by about 10° (from more northerly to northwesterly) with introduction of the IBM slab from 55 to 50 Ma. Most of the required change in absolute Pacific Plate motion must come from elsewhere, and Hu et al. (2022) show that this can be accomplished with termination of intra-oceanic subduction in the north Pacific between 55 and 50 Ma. Other west Pacific subduction zones forming during the Cenozoic, however, may have initially been induced, including the Tonga-Kermadec subduction zone that underwent vertical motions and folding within the future back-arc region prior to and contemporaneously with subduction initiation, indicating a strong compressive force during initiation (Sutherland et al., 2020). Strike–slip motion between Pacific and West Phillipine Sea Plates, along with compression along this boundary due to rearrangement of Pacific Plate driving forces can lead to IBM subduction initiation. The models favor an induced mode of subduction initiation, as the minimum stress (Supplementary Figure S3 solid lines) never vanishes independent of the magnitude of V_SS.

The Vanuatu subduction zone initiated between 12 and 10 Ma, following the collision of the Ontong-Java Plateau with the Vitiaz subduction zone at 16 Ma (Mann and Taira, 2004) and may be an example of induced subduction polarity reversal (Auzende et al., 1988; Yang, 2022). Soon after Vanuatu subduction initiation, the Fiji basin opens at ∼8 Ma, indicating a short initiation time (Auzende et al., 1988; Lallemand and Arcay, 2021). The plate reconstructions show a strong component of strike-slip motion (5–8 cm/yr), Figure 1B dashed blue curve) between the Australian and Pacific plates along the Vanuatu incipient boundary (Figures 1B, D) The strike-slip mechanical model yields a small value for t_SI ($\approx 3$ Myr, Figure 3), because of the large strike-slip velocity V_SS, large convergent velocity V_x and the young plate age (Lallemand and Arcay, 2021).

5 Conclusion

Three Cenozoic subduction zones—Puysegur, IBM and Vanuatu—all show substantial components of strike–slip motion prior to and during initiation. Application of the mechanical model developed here using plate reconstruction constraints all yield estimates of initiation times consistent with values estimated independently from geological observations. Together, the models and observations suggest that subduction initiation can be triggered when margins become progressively weakened to the point that the resisting forces become smaller than the driving forces. Despite the strike-slip velocity having a relatively minor influence on the evolution of plate stress and initiation time compared to the plate convergent velocity, it can still dramatically lower the force required to induce subduction initiation, thereby providing a favorable condition for subduction initiation.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material. The datasets for this study can be found in the CaltechDATA at https://doi.org/110.22002/ae55y-m7q16.

Author contributions

YL contributes to study conception, model setup, manuscript preparation, interpretation MG contributes to study conception, manuscript preparation, interpretation. All authors reviewed the results and approved the final version of the manuscript.

Funding

Supported by the National Science Foundation through awards OCE–1654766, EAR–1645775, and OCE–2049086. This work used Stampede-2 at the Texas Advanced Computer Center from the Extreme Science and Engineering Discovery Environment (XSEDE), which was supported by National Science Foundation.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2023.1156034/full#supplementary-material

References

Arculus, R. J., Gurnis, M., Ishizuka, O., Reagan, M. K., Pearce, J. A., and Sutherland, R. (2019). How to create new subduction zones: A global perspective. Oceanography 32, 160–174. doi:10.5670/oceanog.2019.140

CrossRef Full Text | Google Scholar

Auzende, J.-M., Lafoy, Y., and Marsset, B. (1988). Recent geodynamic evolution of the north Fiji basin (southwest Pacific). Geology 16, 925–929. doi:10.1130/0091-7613(1988)016<0925:rgeotn>2.3.co;2

CrossRef Full Text | Google Scholar

Beavan, J., Samsonov, S., Denys, P., Sutherland, R., Palmer, N., and Denham, M. (2010). Oblique slip on the Puysegur subduction interface in the 2009 july Mw 7.8 Dusky Sound earthquake from GPS and InSAR observations: Implications for the tectonics of southwestern New Zealand. Geophys. J. Int. 183, 1265–1286. doi:10.1111/j.1365-246X.2010.04798.x

CrossRef Full Text | Google Scholar

Bercovici, D., and Ricard, Y. (2014). Plate tectonics, damage and inheritance. Nature 508, 513–516. doi:10.1038/nature13072

PubMed Abstract | CrossRef Full Text | Google Scholar

Bott, M. H. P. (1991). Ridge push and associated plate interior stress in normal and hotspot regions. Tectonophysics 200, 17–32. doi:10.1016/0040-1951(91)90003-b

CrossRef Full Text | Google Scholar

Brune, J. N., Henyey, T. L., and Roy, R. F. (1969). Heat flow, stress, and rate of slip along the San Andreas Fault, California. J. Geophys. Res. 74, 3821–3827. doi:10.1029/JB074i015p03821

CrossRef Full Text | Google Scholar

Choi, H., Kim, S.-S., Dyment, J., Granot, R., Park, S.-H., and Hong, J. K. (2017). The kinematic evolution of the Macquarie plate: A case study for the fragmentation of oceanic lithosphere. Earth Planet. Sci. Lett. 478, 132–142. doi:10.1016/j.epsl.2017.08.035

CrossRef Full Text | Google Scholar

Collot, J. Y., Lamarche, G., Wood, R. A., Delteil, J., Sosson, M., Lebrun, J. F., et al. (1995). Morphostructure of an incipient subduction zone along a transform plate boundary - Puysegur Ridge and Trench. Geology 23, 519–522. doi:10.1130/0091-7613(1995)023<0519:moaisz>2.3.co;2

CrossRef Full Text | Google Scholar

Crameri, F., Magni, V., Domeier, M., Shephard, G. E., Chotalia, K., Cooper, G., et al. (2020). A transdisciplinary and community-driven database to unravel subduction zone initiation. Nat. Commun. 11, 3750–3814. doi:10.1038/s41467-020-17522-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Faulkner, D. R., Mitchell, T. M., Jensen, E., and Cembrano, J. (2011). Scaling of fault damage zones with displacement and the implications for fault growth processes. J. Geophys. Res. Solid Earth 116, B05403. doi:10.1029/2010JB007788

CrossRef Full Text | Google Scholar

Gurnis, M. (2023). “An evolutionary perspective on subduction initiation,” in Dynamics of plate tectonics and mantle convection. Editor J. C. Duarte (Amsterdam: Elsevier), 357–383. doi:10.1016/B978-0-323-85733-8.00003-2

CrossRef Full Text | Google Scholar

Gurnis, M., Hall, C., and Lavier, L. (2004). Evolving force balance during incipient subduction. Geochem. Geophys. Geosystems 5. doi:10.1029/2003gc000681

CrossRef Full Text | Google Scholar

Hall, C., Gurnis, M., Sdrolias, M., Lavier, L. L., and Müller, R. D. (2003). Catastrophic initiation of subduction following forced convergence across fracture zones. Earth Planet. Sci. Lett. 212, 15–30. doi:10.1016/S0012-821X(03)00242-5

CrossRef Full Text | Google Scholar

Hall, R., Ali, J. R., Anderson, C. D., and Baker, S. J. (1995). Origin and motion history of the Philippine Sea Plate. Tectonophysics 251, 229–250. doi:10.1016/0040-1951(95)00038-0

CrossRef Full Text | Google Scholar

Hightower, E., Gurnis, M., and Van Avendonk, H. (2020). A bayesian 3-D linear gravity inversion for complex density distributions: Application to the Puysegur subduction system. Geophys. J. Int. 223, 1899–1918. doi:10.1093/gji/ggaa425

CrossRef Full Text | Google Scholar

Hirth, G., and Kohlstedt, D. L. (1996). Water in the oceanic upper mantle: Implications for rheology, melt extraction and the evolution of the lithosphere. Earth Planet. Sci. Lett. 144, 93–108. doi:10.1016/0012-821x(96)00154-9

CrossRef Full Text | Google Scholar

Hu, J., Gurnis, M., Rudi, J., Stadler, G., and Müller, D. (2022). Dynamics of the abrupt change in Pacific Plate motion around 50 million years ago. Nat. Geosci. 15, 74–78. doi:10.1038/s41561-021-00862-6

CrossRef Full Text | Google Scholar

Hu, J., and Gurnis, M. (2020). Subduction duration and slab dip. Geochem. Geophys. Geosystems 21. doi:10.1029/2019GC008862

CrossRef Full Text | Google Scholar

Humphreys, E. D., and Clayton, R. W. (1990). Tomographic image of the southern California mantle. J. Geophys. Res. 95, 19725. doi:10.1029/jb095ib12p19725

CrossRef Full Text | Google Scholar

Ishizuka, O., Hickey-Vargas, R., Arculus, R. J., Yogodzinski, G., Savov, I. P., Kusano, A., et al. (2018). Age of Izu Bonin mariana arc basement. Earth Planet. Sci. Lett. 481, 80–90. doi:10.1016/j.epsl.2017.10.023

CrossRef Full Text | Google Scholar

Ismail-Zadeh, A., and Tackley, P. J. (2010). Computational methods for geodynamics. Cambridge, U.K.: Cambridge University Press.

Google Scholar

Kohlstedt, D. L., Evans, B., and Mackwell, S. J. (1995). Strength of the lithosphere: Constraints imposed by laboratory experiments. J. Geophys. Res. 100 (17), 17587–17602. doi:10.1029/95jb01460

CrossRef Full Text | Google Scholar

Lallemand, S., and Arcay, D. (2021). Subduction initiation from the earliest stages to self-sustained subduction: Insights from the analysis of 70 Cenozoic sites. Earth-Science Rev. 221, 103779. doi:10.1016/j.earscirev.2021.103779

CrossRef Full Text | Google Scholar

Lamarche, G., and Lebrun, J.-F. (2000). Transition from strike-slip faulting to oblique subduction: Active tectonics at the Puysegur margin, south New Zealand. Tectonophysics 316, 67–89. doi:10.1016/s0040-1951(99)00232-2

CrossRef Full Text | Google Scholar

Lebrun, J.-F., Lamarche, G., and Collot, J.-Y. (2003). Subduction initiation at a strike-slip plate boundary: The Cenozoic Pacific-Australian plate boundary, south of New Zealand. J. Geophys. Res. Solid Earth 108. doi:10.1029/2002jb002041

CrossRef Full Text | Google Scholar

Leng, W., and Gurnis, M. (2011). Dynamics of subduction initiation with different evolutionary pathways. Geochem. Geophys. Geosys. 12. doi:10.1029/2011gc003877

CrossRef Full Text | Google Scholar

Leng, W., and Gurnis, M. (2015). Subduction initiation at relic arcs. Geophys. Res. Lett. 42, 7014–7021. doi:10.1002/2015gl064985

CrossRef Full Text | Google Scholar

Li, H.-Y., Li, X., Ryan, J. G., Zhang, C., and Xu, Y.-G. (2022). Boron isotopes in boninites document rapid changes in slab inputs during subduction initiation. Nat. Commun. 13, 993–1010. doi:10.1038/s41467-022-28637-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Y., and Gurnis, M. (2022). A simple force balance model of subduction initiation. Geophys. J. Int. 232, 128–146. doi:10.1093/gji/ggac332

CrossRef Full Text | Google Scholar

Mann, P., and Taira, A. (2004). Global tectonic significance of the Solomon islands and ontong java plateau convergent zone. Tectonophysics 389, 137–190. doi:10.1016/j.tecto.2003.10.024

CrossRef Full Text | Google Scholar

Mansour, J., Kaluza, O., Giordani, J., Beucher, R., Farringtona, R., Kennedy, G., et al. (2019). underworldcode/underworld2: v2.8.1b. doi:10.5281/zenodo.3384283

CrossRef Full Text | Google Scholar

Maunder, B., Prytulak, J., Goes, S., and Reagan, M. (2020). Rapid subduction initiation and magmatism in the Western Pacific driven by internal vertical forces. Nat. Commun. 11, 1874. doi:10.1038/s41467-020-15737-4

PubMed Abstract | CrossRef Full Text | Google Scholar

McKenzie, D. (1977). “The initiation of trenches: A finite amplitude instability,” in Island arcs, deep sea trenches and back-arc basins, Maurice Ewing Series. Editors M. Talwani, and W. Pitman (Washington, DC: American Geophysical Union), 1, 57–61.

CrossRef Full Text | Google Scholar

Moresi, L., Dufour, F., and Muhlhaus, H. B. (2003). A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. J. Comput. Phys. 184, 476–497. doi:10.1016/s0021-9991(02)00031-1

CrossRef Full Text | Google Scholar

Moresi, L., Gurnis, M., and Zhong, S. (2000). Plate tectonics and convection in the Earth’s mantle: Toward a numerical simulation. Comp. Sci. Engin. 2, 22–33. doi:10.1109/5992.841793

CrossRef Full Text | Google Scholar

Müller, R. D., Zahirovic, S., Williams, S. E., Cannon, J., Seton, M., Bower, D. J., et al. (2019). A global plate model including lithospheric deformation along major rifts and orogens since the Triassic. Tectonics 38, 1884–1907. doi:10.1029/2018tc005462

CrossRef Full Text | Google Scholar

Nikolaeva, K., Gerya, T., and Marques, F. (2010). Subduction initiation at passive margins: Numerical modeling. J. Geophys. Res. 115, B03406. doi:10.1029/2009JB006549

CrossRef Full Text | Google Scholar

Patriat, M., Falloon, T., Danyushevsky, L., Collot, J., Jean, M. M., Hoernle, K., et al. (2019). Subduction initiation terranes exposed at the front of a 2 Ma volcanically–active subduction zone. Earth Planet. Sci. Lett. 508, 30–40. doi:10.1016/j.epsl.2018.12.011

CrossRef Full Text | Google Scholar

Qing, J., Liao, J., Li, L., and Gao, R. (2021). Dynamic evolution of induced subduction through the inversion of spreading ridges. J. Geophys. Res. Solid Earth 126, e2020JB020965. doi:10.1029/2020jb020965

CrossRef Full Text | Google Scholar

Reagan, M. K., Heaton, D. E., Schmitz, M. S., Pearce, J. A., Shervais, J. W., and Koppers, A. P. (2019). Forearc ages reveal extensive short–lived and rapid seafloor spreading following subduction initiation. Earth Planet. Sci. Lett. 506, 520–529. doi:10.1016/j.epsl.2018.11.020

CrossRef Full Text | Google Scholar

Savage, H. M., and Brodsky, E. E. (2011). Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones. J. Geophys. Res. Solid Earth 116, B03405. doi:10.1029/2010JB007665

CrossRef Full Text | Google Scholar

Shuck, B., Gulick, S. P. S., Van Avendonk, H. J. A., Gurnis, M., Sutherland, R., Stock, J., et al. (2022). Stress transition from horizontal to vertical forces during subduction initiation. Nat. Geosci. 15, 149–155. doi:10.1038/s41561-021-00880-4

CrossRef Full Text | Google Scholar

Shuck, B., Van Avendonk, H., Gulick, S. P. S., Gurnis, M., Sutherland, R., Stock, J., et al. (2021). Strike-slip enables subduction initiation beneath a failed rift: New seismic constraints from Puysegur margin, New Zealand. Tectonics 40, e2020TC006436. doi:10.1029/2020tc006436

CrossRef Full Text | Google Scholar

Stern, R. J., and Bloomer, S. H. (1992). Subduction zone infancy - examples from the Eocene izu–bonin–mariana and jurassic California arcs. Geol. Soc. Am. Bull. 104, 1621–1636. doi:10.1130/0016-7606(1992)104<1621:szieft>2.3.co;2

CrossRef Full Text | Google Scholar

Sutherland, R., Barnes, P., and Uruski, C. (2006). Miocene-Recent deformation, surface elevation, and volcanic intrusion of the overriding plate during subduction initiation, offshore southern Fiordland, Puysegur margin, southwest New Zealand. N.Z. J. Geol. geophys. 49, 131–149. doi:10.1080/00288306.2006.9515154

CrossRef Full Text | Google Scholar

Sutherland, R., Dickens, G. R., Blum, P., Agnini, C., Alegret, L., Asatryan, G., et al. (2020). Continental-scale geographic change across Zealandia during Paleogene subduction initiation. Geology 48, 419–424. doi:10.1130/G47008.1

CrossRef Full Text | Google Scholar

Thielmann, M., and Kaus, B. J. P. (2012). Shear heating induced lithospheric-scale localization: Does it result in subduction? Earth Planet. Sci. Lett. 359, 1–13. doi:10.1016/j.epsl.2012.10.002

CrossRef Full Text | Google Scholar

Torsvik, T. H., Doubrovine, P. V., Steinberger, B., Gaina, C., and Spakman, W. (2017). Pacific plate motion change caused the Hawaiian-Emperor bend. Nat. Commun. 8, 15660. doi:10.1038/ncomms15660

PubMed Abstract | CrossRef Full Text | Google Scholar

Toth, J., and Gurnis, M. (1998). Dynamics of subduction initiation at preexisting fault zones. J. Geophys. Research-Solid Earth 103, 18053–18067. doi:10.1029/98jb01076

CrossRef Full Text | Google Scholar

Whittaker, J. M., Muller, R. D., Leitchenkov, G., Stagg, H., Sdrolias, M., Gaina, C., et al. (2007). Major Australian–Antarctic plate reorganization at Hawaiian–Emperor bend time. Science 318, 83–86. doi:10.1126/science.1143769

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, G. (2022). Subduction initiation triggered by collision: A review based on examples and models. Earth-Science Rev. 232, 104129. doi:10.1016/j.earscirev.2022.104129

CrossRef Full Text | Google Scholar

Zhong, X., and Li, Z.-H. (2023). Compression at strike-slip fault is a favorable condition for subduction initiation. Geophys. Res. Lett. 50, e2022GL102171. doi:10.1029/2022gl102171

CrossRef Full Text | Google Scholar

Zhong, X., and Li, Z.-H. (2019). Forced subduction initiation at passive continental margins: Velocity-driven versus stress-driven. Geophys. Res. Lett. 46, 11054–11064. doi:10.1029/2019gl084022

CrossRef Full Text | Google Scholar

Zhou, X., Li, Z.-H., Gerya, T. V., Xu, Z., and Zhang, J. (2018). Subduction initiation dynamics along a transform fault control trench curvature and ophiolite ages. Geology 46, 607–610. doi:10.1130/G40154.1

CrossRef Full Text | Google Scholar

Zoback, M. D., Zoback, M. L., Mount, V. S., Suppe, J., Eaton, J. P., Healy, J. H., et al. (1987). New evidence on the state of stress of the San Andreas Fault system. Science 238, 1105–1111. doi:10.1126/science.238.4830.1105

PubMed Abstract | CrossRef Full Text | Google Scholar