- Seismological Laboratory, California Institute of Technology, Pasadena, CA, United States
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, tSI, marking the moment when intraplate stress flips sign, is validated with a computational models. The analytical solution shows that tSI 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?
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1156034/feart-11-1156034-HTML/image_m/feart-11-1156034-g001.jpg)
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 (
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,
Viscosity follows the non-Newtonian Arrhenius law,
where μ is the shear modulus, η the viscosity, T the temperature, and
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.
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.
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, Fxx, and strike-slip velocity, VSS. 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, Fxx, initiation time, tSI, and total work, WSI, 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.,
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1156034/feart-11-1156034-HTML/image_m/feart-11-1156034-g002.jpg)
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 Fxx as a function of plate motion ℓ for models with different strike slip velocity VSS 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 WSI of every numerical model (color-coded solid circle) in the parametric grid search varying the strike slip velocity VSS and ɛP0. Contours show the prediction of WSI from analytical model.
In 2D cross-sectional models (with no strike-slip motion), Fxx can be formulated analytically via a force balance analysis on the subducting plate, such that Fxx, 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
Finally, we evaluate tSI from the analytical model varying Vx, VSS, ɛ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 Vx has a first-order control on the predicted tSI. 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 tSI is dominated by slab pull when plate is young (
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1156034/feart-11-1156034-HTML/image_m/feart-11-1156034-g003.jpg)
FIGURE 3. Analytical prediction of initiation time tSI under different convergent velocity Vx and subducting plate age with (A). VSS = 0, ɛP0 = 1, (B). VSS = 3 cm/yr, (C). VSS = 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 VSS 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 Fxx 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 Fxx. When Fxx 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, VSS, fxx) = (2, 0 cm/yr, 8.57 × 1012 N/m), as a reference (Supplementary Figure S2B) In this case, due to the insufficient compression to induce subduction initiation (Fxx = 8.57 × 1012 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 × 1012 N/m (Bott, 1991; Toth and Gurnis, 1998) and so the Fxx 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 Vss 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.
![www.frontiersin.org](https://www.frontiersin.org/files/Articles/1156034/feart-11-1156034-HTML/image_m/feart-11-1156034-g004.jpg)
FIGURE 4. Fault evolution of a high resolution models with (ɛP0, VSS, Vp) = (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
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 Vx, VSS, and plate age for Puysegur, we predict a range of tSI from the analytical model (Figure 3, blue box). The northern Puysegur trench corresponds to the top right corner of the blue box (Vx = 1 cm/yr, plate age = 25 Ma), yielding a model tSI ≈ 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 tSI (
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 VSS and Vx 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 tSI range from 5 to 10 Myr with VSS =3–6 cm/yr. At 50 Ma, the reconstruction gives a larger range of VSS (0–6 cm/yr), and the analytical model predicts a tSI ≈ 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 40Ar/39Ar 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 VSS.
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 tSI (
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
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
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
Bercovici, D., and Ricard, Y. (2014). Plate tectonics, damage and inheritance. Nature 508, 513–516. doi:10.1038/nature13072
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
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
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
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
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
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
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
Gurnis, M., Hall, C., and Lavier, L. (2004). Evolving force balance during incipient subduction. Geochem. Geophys. Geosystems 5. doi:10.1029/2003gc000681
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
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
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
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
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
Hu, J., and Gurnis, M. (2020). Subduction duration and slab dip. Geochem. Geophys. Geosystems 21. doi:10.1029/2019GC008862
Humphreys, E. D., and Clayton, R. W. (1990). Tomographic image of the southern California mantle. J. Geophys. Res. 95, 19725. doi:10.1029/jb095ib12p19725
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
Ismail-Zadeh, A., and Tackley, P. J. (2010). Computational methods for geodynamics. Cambridge, U.K.: Cambridge University Press.
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
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
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
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
Leng, W., and Gurnis, M. (2011). Dynamics of subduction initiation with different evolutionary pathways. Geochem. Geophys. Geosys. 12. doi:10.1029/2011gc003877
Leng, W., and Gurnis, M. (2015). Subduction initiation at relic arcs. Geophys. Res. Lett. 42, 7014–7021. doi:10.1002/2015gl064985
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
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
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
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
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
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.
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
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
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
Nikolaeva, K., Gerya, T., and Marques, F. (2010). Subduction initiation at passive margins: Numerical modeling. J. Geophys. Res. 115, B03406. doi:10.1029/2009JB006549
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
Keywords: geodynamics, subduction, subduction initiation, faults, plate kinematics
Citation: Li Y and Gurnis M (2023) Strike slip motion and the triggering of subduction initiation. Front. Earth Sci. 11:1156034. doi: 10.3389/feart.2023.1156034
Received: 01 February 2023; Accepted: 30 March 2023;
Published: 17 April 2023.
Edited by:
Taras Gerya, ETH Zürich, SwitzerlandReviewed by:
Gaoxue Yang, Chang’an University, ChinaZhong-Hai Li, University of Chinese Academy of Sciences, China
Jie Liao, School of Earth Sciences and Engineering, Sun Yat-Sen University, China
Copyright © 2023 Li and Gurnis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yida Li, yidali@caltech.edu