Can subduction initiation at a transform fault be spontaneous?

We present an extensive parametric exploration of the feasibility of “spontaneous” subduction initiation, i.e., lithospheric gravitational collapse without any external forcing, at a transform fault (TF). We first seek candidates from recent subduction initiation events at an oceanic TF that could fulfill the criteria of spontaneous subduction and retain three natural cases: Izu–Bonin–Mariana, Yap, and Matthew and Hunter. We next perform an extensive exploration of conditions allowing for the spontaneous gravitational sinking of the older oceanic plate at a TF using 2-D thermomechanical simulations. Our parametric study aims at better delimiting the ranges of mechanical properties necessary to achieve the old plate sinking (OPS). The explored parameter set includes the following: crust and TF densities, brittle and ductile rheologies, and the width of the weakened region around the TF. We focus on characterizing the OPS conditions in terms of (1) the reasonable vs. unrealistic values of the mechanical parameters and (2) a comparison to modern cases of subduction initiation in a TF setting. When modeled, OPS initiates following one of two distinct modes, depending mainly on the thickness of the overlying younger plate. The asthenosphere may rise up to the surface above the sinking old plate, provided that the younger plate remains motionless (verified for ages ≥5 Myr, mode 1). For lower younger plate ages (typically ≤2 Myr), the younger plate is dragged toward the older plate, resulting in a double-sided subduction (mode 2). When triggered, spontaneous OPS is extremely fast. The parameters that exert the strongest control over whether OPS can occur or not are the brittle properties of the shallow part of the lithosphere, which affect the plate resistance to bending, the distance away from the TF over which weakening is expected, and the crust density. We find that at least one mechanical parameter has to be assigned an unrealistic value and at least two other ones must be set to extreme ranges to achieve OPS, which we do not consider realistic. Furthermore, we point out inconsistencies between the processes and consequences of lithospheric instability, as modeled in our experiments and geological observations of subduction infancy, for the three natural candidates of subduction initiation by spontaneous OPS. We conclude that spontaneous instability of the thick older plate at a TF evolving into mature subduction is an unlikely process of subduction initiation in modern Earth conditions.

1 Introduction

The process of spontaneous subduction deserves to be explored again following recent discoveries during Ocean Drilling Project 375 in the backarc of the proto-Izu–Bonin subduction zone. The nature and age of the basaltic crust drilled there appeared to be similar to those of the forearc basalts underlying the boninites of the present Izu–Bonin forearc (Hickey-Vargas et al., 2018). The consequences of this discovery are controversial since they are supposed to support the concept of spontaneous subduction for some authors (Arculus et al., 2015; Stern and Gerya, 2018), whereas for other authors, they do not (Keenan and Encarnación, 2016; Lallemand, 2016).

The notion of “spontaneous subduction” originates from two observations: (1) Uyeda and Kanamori (1979) first described the Mariana-type extreme subduction mode, where an old oceanic plate sunk, driven by its weight excess into a vertical slab in association with backarc extension. (2) A few years later, in the early 1980s, Bonin Island volcanic rock analysis and deep sea drilling (Leg 60) in the adjacent Izu–Bonin–Mariana (IBM) subduction zone forearc revealed rocks called boninites that combined the characteristics of arc-lavas and MORB (Mid-ocean ridge basalt; Natland and Tarney, 1981; Bloomer and Hawkins, 1983). A conceptual model was then proposed by Stern and Bloomer (1992) reconciling these observations, in which an old plate may sink in the mantle under its weight along the weak boundary formed by a transform fault. Numerical models (Hall and Gurnis, 2003; Gurnis et al., 2004) first failed to support this process of spontaneous subduction and concluded that a tectonic force was required to initiate subduction. Later, they finally succeeded in simulating spontaneous subduction in specific contexts, such as lithospheric collapse around a plume head (Whattam and Stern, 2015) or the conjunction of a large density contrast with a very weak fault zone between the adjacent lithospheres (Leng and Gurnis, 2015).

In this study, we will adopt the definition of Stern and Gerya (2018): spontaneous subduction is caused by forces originating at the subduction initiation site and not elsewhere (Fig. 1b). They define three different settings where spontaneous subduction may develop: passive margin, transform fault (TF), or plume head. The only Cenozoic examples that were attributed by Stern and Gerya (2018) to potential sites of spontaneous subduction initiation, i.e., IBM and Tonga-Kermadec, correspond to the TF setting (Fig. 1a). In these two examples, the relics of the subduction initiation stage date back to the Eocene and are thus subject to controversy. We first recall the natural examples for which oceanic TFs or fracture zones might have evolved into a subduction zone. Then, numerical models addressing subduction initiation processes in a similar context are analyzed before developing our own numerical approach. The range of parameters allowing for spontaneous subduction initiation in our models will finally be compared with the reasonable values characterizing the natural processes.

Table 1 Oceanic TFs or fracture zones, where potential subduction initiated during Cenozoic times.

References: 1: Stern and Bloomer (1992); 2: Lallemand (2016); 3: Sdrolias and Müller (2006); 4: Gaina and Müller (2007); 5: Fujiwara et al. (2000); 6: Hawkins and Batiza (1977); 7: Hegarty and Weissel (1988); 8: Patriat et al. (2015); 9: Mortimer et al. (2014); 10: Hegarty et al. (1983); 11: Lebrun et al. (2003); 12: Meckel et al. (2005); 13: Bonatti et al. (1996); 14: Meckel et al. (2005); 15: Deschamps and Lallemand (2002); 16: Abecassis et al. (2016); 17: Deschamps et al. (2000); 18: Sibuet et al. (2002); 19: Hilde and Lee (1984); 20: Doo et al. (2015); 21: Kuo et al. (2009); 22: Eakin et al. (2015); 23: Patriat et al. (2011); 24: Pichot et al. (2012); 25: Gutscher et al. (2002); 26: Duarte et al. (2013); 27: Tortella et al. (1997); 28: Maia et al. (2016); 29: Fournier et al. (2011).

Figure 1 Various tectonic settings leading to vertical motion and/or convergence at transform plate boundaries, as detailed in Table 1. The convergent black heavy arrows represent far-field tectonic forces. The red light arrows outline the sense of motion of one plate with respect to the other. The red crosses and dots in circles indicate transform motion. The thicker plate is the older one.

1.1 From oceanic transform faults or fracture zones to subduction in nature

Table 1 and Fig. 1 summarize Cenozoic settings where oceanic TFs or fracture zones underwent deformation that sometimes evolved into subduction and at other times did not. The regions are classified in Table 1 such that the older plate (OP) underthrusts the younger in the first group (Fig. 1b, c, d, IBM, Yap, Matthew and Hunter, Mussau, Macquarie, and Romanche), the downgoing plate is the youngest in the second group (Fig. 1d, e, Hjort, Gagua, Barracuda and Tiburon), and finally those for which it appears to be impossible to determine the relative age of one plate with respect to the other at the time of initiation (Fig. 1f, Gorringe, St Paul and Owen). The analysis of all these natural cases shows that the 3-D setting and far-field boundary conditions are likely to play a major role in subduction initiation and on the selected age (old/young) of the subducting plate. Earlier studies showed that compression prevailed in the upper plate at the time of initiation for most of them, while it is unknown for IBM and Yap. In these two regions, subduction started more than 20 Myr ago (Hegarty and Weissel, 1988; Ishizuka et al., 2011), but, soon after they were initiated, they underwent one of the strongest episodes of subduction erosion on Earth (Natland and Tarney, 1981; Hussong and Uyeda, 1981; Bloomer, 1983; Lallemand, 1995), so all remnants of their forearc at the time of initiation were consumed (Lallemand, 2016, and references therein). Geological evidence of the stress state at initiation is thus either subducted or deeply buried beneath the remnant Palau-Kyushu Ridge. To date, some authors (e.g., Ishizuka et al., 2018; Stern and Gerya, 2018) still argue that spreading, i.e., extension, occurred over a broad area from the backarc to the forearc at the time of subduction initiation. Backarc extension concomitant with subduction initiation under compressive stress is compatible, as exemplified by the recent case of Matthew and Hunter at the southern termination of the New Hebrides subduction zone (Patriat et al., 2015, Fig. 1d). There, the authors suggest that the collision of the Loyalty Ridge with the New Hebrides Arc induced the fragmentation of the North Fiji Basin (Eissen spreading center and Monzier rift), whose extension yielded, in turn, a compressive stress along the southern end of the transform boundary (or STEP fault), accommodating the trench rollback of the New Hebrides trench. It is important to note that the geodynamic context of the Matthew and Hunter region is very similar to the one of the IBM protosubduction (Deschamps and Lallemand, 2002, 2003; Patriat et al., 2015; Lallemand, 2016). Rifting and spreading in a direction normal to the TF has been documented at the time of subduction initiation. Since the conditions of spontaneous subduction do not require compressive stress, but rather the sinking of the oldest plate under its weight excess, and because of the lack of geological records of what happened there, we consider that IBM and Yap subduction initiation might be either spontaneous (Fig. 1b) or forced (Fig. 1c). To decipher between these two hypotheses, we conduct a series of numerical simulations.

1.2 Modeling of spontaneous subduction initiation at a transform fault in previous studies

Numerical experiments have shown that old plate sinking (OPS) could spontaneously occur for a limited viscosity contrast between lithospheres and the underlying asthenosphere (Matsumoto and Tomoda, 1983) in a model neglecting thermal effects. However, without imposed convergence, subduction initiation failed when thermal diffusion was taken into account, even in the most favorable case of an old and thick plate facing a section of asthenosphere (Hall and Gurnis, 2003; Baes and Sobolev, 2017), unless the density offset at the TF was emphasized by including a thick and buoyant crust at the younger plate (YP) surface (Leng and Gurnis, 2015). In most cases showing the instability of the thick plate, lateral density contrasts at the TF are maximized by imposing at the TF an extremely thin younger plate (0 or 1 Myr old at the location where instability initiates) in front of a thicker plate, whose age is chosen between 40 and 100 Myr, either in 2-D (Nikolaeva et al., 2008) or 3-D (Zhu et al., 2009, 2011; Zhou et al., 2018). For similar plate age pairs, Gerya et al. (2008) showed that successful spontaneous initiation requires the OP slab surface to be sufficiently lubricated and strongly weakened by metasomatism to decouple the two adjacent plates as plate sinking proceeds, while the dry mantle is supposed to be moderately resistant to bending. Assuming such “weak” rheological structure, OPS triggering occurs and results in an asthenosphere rise in the vicinity of the subduction hinge, which yields a fast spreading (from a few centimeters per year to >1 m yr⁻¹). It has been described as a “catastrophic” subduction initiation (Hall and Gurnis, 2003). This catastrophic aspect is hampered when thicker YPs are considered (10 to 20 Myr old), when crustal and mantle rheologies are less weak, and when shallow plate weakening develops progressively through time, e.g., by pore fluid pressure increase with sea water downward percolation in a low-permeability matrix (Dymkova and Gerya, 2013).

These previous numerical studies have helped to unravel the conditions leading to OPS without any imposed external forcing. Nevertheless, recent incipient subduction zones, the most likely to correspond to initiation by spontaneous sinking at a TF, are not all associated with a significant plate age offset at plate boundaries (Matthew and Hunter, Yap, Table 1). We thus propose a new investigation of the conditions of OPS to address the following three questions. What are the mechanical parameter ranges allowing for OPS, especially for the TF settings that are the closest to spontaneous subduction conditions? Are these parameter ranges reasonable? Are the modeled kinetics and early deformation compatible with natural cases observations?

We choose a simplified setup, without fluid percolation simulations and in 2-D, to allow for a broad parameter exploration with an accurate numerical resolution.

Table 2 Constant names and values.

2 Model setup

The numerical model solves the momentum, energy, and mass conservation equations, assuming that rocks are incompressible, except for the thermal buoyancy term in the momentum equation and for the adiabatic heating term in the energy equation (extended Boussinesq approximation). As shear heating has been shown to significantly improve strain localization within the subduction interface (Doin and Henry, 2001; Thielmann and Kaus, 2012), it is included in the heat conservation equation, as well as a uniform heat production (Table 2). The simulation box, 2220 km wide and 555 km thick, is chosen to be large enough to simulate convective interactions between shallow lithospheres and deep mantle that may be involved in the process of subduction initiation (Fig. 2).

Density (ρ) is assumed to be temperature and composition dependent:

$$\begin{array}{}\text{(1)}& \mathit{\rho }(C,T)={\mathit{\rho }}^{\mathrm{ref}}\left(C\right)(\mathrm{1}-\mathit{\alpha }(T-T_{\mathrm{s}}\left)\right),\end{array}$$

where ρ^ref is the reference density at the surface, C is composition (mantle, oceanic crust, or weak material; Sect. 2.3), α is the thermal expansion coefficient, T is temperature, and T_s is the surface temperature (Table 2). For the mantle, ${\mathit{\rho }}_{\mathrm{m}}^{\mathrm{ref}}$ is fixed to 3300 kg m⁻³, while ρ^ref for the oceanic crust and the weak material is varied from one experiment to another (Sect. 2.4).

Figure 2 (a) Boundary conditions, initial thermal state, and material distribution at simulation start. One white isotherm every 200 ^∘C. If the box bottom is open, a vertical resistance against flow is imposed to simulate a viscosity jump approximately 10 times higher below the simulation box (Ribe and Christensen, 1994; Arcay, 2012). The red dotted line represents the hot thermal anomaly ($\mathrm{\Delta }T=+\mathrm{250}$ ^∘C) imposed in some experiments. (b) Close-up on the TF structure. L_w is the width at the surface of the younger plate and of the older plate (aged of A_y and of A_o Myr, respectively) over which the oceanic crust is assumed to have been altered and weakened by the TF activity. The meaning of labels 1 to 4 is given in Sect. 2.3.

2.1 Rheology

We combine a pseudo-brittle rheology to a non-Newtonian ductile law. Pseudo-brittle rheology is modeled using a yield stress, τ_y, increasing with depth, z:

$$\begin{array}{}\text{(2)}& {\mathit{\tau }}_{\mathrm{y}}=C_{\mathrm{0}}+\mathit{\gamma }\left(C\right)\mathit{\rho }gz,\end{array}$$

where C₀ is the cohesive strength at the surface (Table 2), γ is a function of composition C, ρ is density, and g is the gravity acceleration. The parameter γ represents the yield strength increase with depth and can be related to the coefficient of internal friction of the Coulomb–Navier criterion (Sect. 2.5). To simplify, we tag γ as the brittle parameter. The relationship between the lithostatic pressure ρgz and the normal stress σ_n applied on the brittle fault will be derived in Sect. 2.5.1. The brittle deviatoric strain rate is computed assuming the following relationship (Doin and Henry, 2001): $\dot{\mathit{\epsilon }}={\dot{\mathit{\epsilon }}}_{\mathrm{ref}}(\mathit{\tau }/{\mathit{\tau }}_{\mathrm{y}}{)}^{n_p}$, where $\dot{\mathit{\epsilon }}$ is the second invariant of the deviatoric strain rate tensor, ${\dot{\mathit{\epsilon }}}_{\mathrm{ref}}$ is a reference strain rate, and n_p is a large exponent (Table 2). In the plastic domain, strain rates are close to zero if τ≪τ_y but become very large as soon as stress exceeds the yield stress τ_y. Recalling that $\mathit{\tau }=\mathit{\nu }\dot{\mathit{\epsilon }}$, the plastic viscosity, ν_b, is written as follows:

$$\begin{array}{}\text{(3)}& {\mathit{\nu }}_b={\mathit{\tau }}_{\mathrm{y}}{\displaystyle \frac{{\dot{\mathit{\epsilon }}}^{\frac{\mathrm{1}}{n_p}-\mathrm{1}}}{{\dot{\mathit{\epsilon }}}_{\mathrm{ref}}^{\frac{\mathrm{1}}{n_p}}}}.\end{array}$$

A dislocation creep rheology is simulated using a non-Newtonian viscosity ν_d, defined by

$$\begin{array}{}\text{(4)}& {\mathit{\nu }}_{\mathrm{d}}=B_{\mathrm{0}}\exp \left({\displaystyle \frac{E_{\mathrm{a}}\left(C\right)+V_{\mathrm{a}}\mathit{\rho }gz}{nRT}}\right){\dot{\mathit{\epsilon }}}^{\frac{\mathrm{1}}{n}-\mathrm{1}},\end{array}$$

where B₀ is a pre-exponential factor, E_a is the activation energy depending on composition C, V_a is the activation volume, n is the non-Newtonian exponent, and R is the ideal gas constant (Table 2). The effective viscosity ν_eff is computed assuming that the total deformation is the sum of brittle and ductile deformations. Note that the brittle behavior acts as a maximum viscosity cutoff. Regarding strain rate, a minimum cutoff is set to $\mathrm{2.6}\times {\mathrm{10}}^{-\mathrm{21}}$ s⁻¹, but no maximum cutoff is imposed.

Table 3 List of simulations quoted in the text. See the data in the Supplement for the complete simulation list.

^a If one value only is indicated, the oceanic crust (3) in Fig. 2 is assumed to have the same brittle parameter as the weak material forming domains (1) and (2). ^b If only one value indicated, then ρ_c=ρ_TF. ^c When one value only is indicated, L_w is identical on both plates. ^d The weak material is imposed at the fault zone (1) only (Fig. 2). “tw”: thermal transition width at the plate boundary. ΔT_p: temperature anomaly within the plume head. ν_ref: reference viscosity at the lithosphere–asthenosphere boundary (2.74×10¹⁹ Pa s⁻¹). OPS: older plate sinking. YPVSI: YP vertical subduction initiation (as in Fig. 5.6). YP retreat: backward drift of the younger plate, as sketched in Fig. 4.3. YPS: younger plate sinking, as sketched in Fig. 4.4b. Double SI: double subduction initiation (as in Fig. 4.5). SB: slab break off. YPD: young plate dragging and sinking into the mantle.

2.2 Initial thermal structure and boundary conditions

We investigate a wide range of lithosphere age pairs, the younger plate (YP) age, A_y, varying from 0 to 40 Myr, and the older plate (OP) age, A_o, from 5 to 150 Myr (Table 3), to cover the plate age ranges observed in nature (Table 1). The thickness of a lithosphere is here defined by the depth of the 1200 ^∘C isotherm, z_LB(A), classically estimated using the half-space cooling model (Turcotte and Schubert, 1982) by

$$\begin{array}{}\text{(5)}& z_{\mathrm{LB}}\left(A\right)=\mathrm{2}{\mathrm{erf}}^{-\mathrm{1}}\left(\mathrm{0.9}\right)\sqrt{\mathit{\kappa }A}\sim \mathrm{2.32}\sqrt{\mathit{\kappa }A},\end{array}$$

where κ is the thermal diffusivity (Table 2) and A is the plate age. However, the half-space cooling model, as well as some variations of it such as the global median heat flow model (GDH1; Stein and Stein, 1992), have been questioned (Doin et al., 1996; Dumoulin et al., 2001; Hasterok, 2013; Qiuming, 2016). Indeed, such conductive cooling models predict too cold young oceanic plates (by ∼100 to 200 ^∘C) compared to the thermal structure inferred from high-resolution shear wave velocities, such as in the vicinity of the East Pacific Rise (Harmon et al., 2009). Similarly, worldwide subsidence of young seafloors is best modeled by taking into account, in addition to a purely lithosphere conductive cooling model, a dynamic component, likely related to the underlying mantle dynamics (Adam et al., 2015). Recently, Grose and Afonso (2013) have proposed an original and comprehensive model for oceanic plate cooling, which accurately reproduces the distribution of heat flow and topography as a function of seafloor age. This approach leads to young plates (<50 Myr) 100 to 200 ^∘C hotter than predicted using the half-space cooling and Parsons and Sclater (1977) models, especially in the shallowest part of the lithosphere. This discrepancy notably comes from, first, heat removal in the vicinity of the ridge by hydrothermal circulation and, second, the presence of an oceanic crust on top of the lithospheric mantle that insulates it from the cold (0 ^∘C) surface and slows down its cooling and thickening. Taking into account these two processes reduces the surface heat flows predicted by the GDH1 model by 75 % (Grose and Afonso, 2013). Our study focuses on young oceanic plates that are the most frequent at TFs (A_y≲60 Myr, Table 1). Therefore, we calculate lithospheric thicknesses z_LB(A) as 0.75 of the ones predicted by half-space cooling model. Moreover, plates warmer than predicted by the half-space cooling model are consistent with the hypothesis of small-scale convection occurring at the base of very young oceanic lithospheres, i.e., younger than a threshold encompassed between 5 and 35 Myr (Buck and Parmentier, 1986; Morency et al., 2005; Afonso et al., 2008). An early small-scale convection process would explain short-wavelength gravimetric undulations in the plate motion direction in the central Pacific and east-central Indian oceans detected at plate ages older than 10 Myr (e.g., Haxby and Weissel, 1986; Cazenave et al., 1987). Buck and Parmentier (1986) have shown that the factor ${\mathrm{erf}}^{-\mathrm{1}}\left(\mathrm{0.9}\right)\sim \mathrm{1.16}$ in Eq. (5) must be replaced by a value encompassed between 0.74 and 0.93 to fit the plate thicknesses simulated when early small-scale convection is modeled, depending on the assumed asthenospheric viscosity. This is equivalent to applying a corrective factor between 0.74/1.16–0.64 and 0.93/1.16–0.80, which is consistent with the lithosphere thicknesses inferred from heat flow modeling by Grose and Afonso (2013). Between the surface and z_LB(A), the thermal gradient is constant.

The transform fault, located at the middle of the box top (x=1110 km), is modeled by a stair-step function joining the isotherms of the adjacent lithospheres (Fig. 2). We test the effect of the TF thermal state, which should be cooled by conduction in the case of an inactive fracture zone, in a few simulations (Sect. 3.3).

Moreover, we test the possible influence of the asthenospheric thermal state at initiation, either uniform over the whole box or locally marked by thermal anomalies resulting from the small-scale convection observed in a preliminary computation of mantle thermal equilibrium (Fig. 2). The results show that the process of subduction initiation, in the case of success or failure, does not significantly depend on the average asthenospheric thermal structure. Nevertheless, in a few experiments, we impose at the start of simulation a thermal anomaly mimicking a small plume head ascending right below the TF, 200 km wide and ∼75 km high, whose top is located at 110 km depth at start of simulation (Fig. 2). The plume thermal anomaly ΔT_plume is set to 250 ^∘C (Table 3). Regarding boundary conditions, slip is free at the surface and along vertical sides. We test the effect of the box bottom condition, either closed and free-slip or open to mantle in- and outflows. When the box bottom is open, a vertical resistance against flow is imposed along the box base, mimicking a viscosity jump 10 times higher than above (Ribe and Christensen, 1994; Arcay, 2017). The results show that the bottom mechanical condition does not modify the future evolution of the fracture zone. The thermal boundary conditions are depicted in Fig. 2.

2.3 Lithological structure at simulation start

The TF lithological structure is here simplified by considering three different lithologies only: the vertical layer forming the fault zone between the two oceanic lithospheres (label 1 in Fig. 2) and assumed to be the weakest material in the box, the oceanic crust (label 3), and the mantle (label 4). In all experiments, the Moho depth is set to 8.3 km for both oceanic lithospheres, and the width of the vertical weak zone forming the fault 1 is equal to 8.3 km. The depth of the weak vertical zone 1 depends on the chosen older plate age, A_o: it is adjusted to be a bit shallower than the OP base, by ∼15 to 30 km. Furthermore, we want to test the effect of the lateral extent of this weakening, outside the gouge fault, L_w (label 2 in Fig. 2). Indeed, depending on the type of TF, the weak zone width may be limited to ∼8 km, such as for the Discovery and Kane faults (Searle, 1983; Detrick and Purdy, 1980; Wolfson-Schwehr et al., 2014), implying L_w=0 km in our model or, in contrast, that the weak zone width may reach 20 to 30 km, such as for the Quebrada or Gofar TFs (Searle, 1983; Fox and Gallo, 1983); thus, L_w can be varied up to 22 km. In most experiments, we impose the same value for the lateral extent of crust weakening on both lithospheres: L_w(A_o)=L_w(A_y), except in a few simulations.

2.4 Parametric study derived from force balance

The first-order forces driving and resisting subduction initiation at a transform fault indicate which mechanical parameters would be worth testing to study OPS triggering. Without any external forcing, the unique driving force to consider is (1) the plate weight excess relative to the underlying mantle. Subduction is hampered by (2) plate resistance to deformation and bending; (3) the TF resistance to shearing; and (4) the asthenosphere strength, resisting plate sinking (e.g., McKenzie, 1977; Cloetingh et al., 1989; Mueller and Phillips, 1991; Gurnis et al., 2004). We vary the mechanical properties of the different lithologies forming the TF area to alter the incipient subduction force balance. The negative plate buoyancy (1) is related to the plate density, here dependent only on the thermal structure and plate age A (Sect. 2.2) since we do not explicitly model density increase of metamorphized (eclogitized) oceanic crust. Nonetheless, we vary the crust density, ρ_c, imposed at the start of simulation along the plate surface to test the potential effect on plate sinking. We also investigate how the density of the weak layer forming the interplate contact, ρ_TF, which is not well known, may either resist plate sinking (if buoyant) or promote it (if dense). The plate strength and flexural rigidity (2) are varied in our model by playing on different parameters. First, we test the rheological properties of the crustal layer both in the brittle and ductile realms, by varying γ_c and $E_{\mathrm{a}}^{\mathrm{c}}$ (Eqs. 2 and 4). Second, the lithospheric mantle strength is varied through the mantle brittle parameter, γ_m, that controls the maximum lithospheric stress in our model. Third, we vary the lateral extent (L_w) of the shallow lithosphere weakened domain, related to the crust alteration likely to occur in the vicinity of the TF.

We study separately the influence of these six mechanical parameters (ρ_c, ρ_TF, γ_c, $E_{\mathrm{a}}^{\mathrm{c}}$, γ_m, L_w) for most plate age pairs. The TF strength (3) is often assumed to be quite low at the interplate contact (Gurnis et al., 2004; Gerya et al., 2008). We thus fill the TF “gouge” with the weak material (labeled 1 in Fig. 2) and, in most experiments, set it as ${\mathit{\gamma }}_{\mathrm{TF}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$. In some experiments, we replace the weak material filling the TF gouge by the more classical oceanic crust (labeled 3 in Fig. 2) to test the effect of a stiffer fault. In that case, ${\mathit{\gamma }}_{\mathrm{TF}}={\mathit{\gamma }}_{\mathrm{c}}=\mathrm{0.05}$ and L_w=0 km: the TF and both plate surfaces are made of gabbroic oceanic crust (Table 3). Note that when ${\mathit{\gamma }}_{\mathrm{c}}={\mathit{\gamma }}_{\mathrm{TF}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$, the weak layer and the oceanic crust are mechanically identical, and the weak layer then entirely covers the whole plate surface (L_w=1100 km). Similarly, as the activation energy $E_{\mathrm{a}}^{\mathrm{c}}$ is the same for the oceanic crust and the weak material, assuming a low ductile strength for the TF is equivalent to covering the whole plate surface by the weak layer (setting L_w=1100 km).

Apart from the six main physical properties that are repeatedly tested (Sect. 2.5), we perform additional experiments for a limited number of plate age combinations to investigate a few extra parameters. In this set of simulations, we vary the asthenosphere resistance competing against plate sinking (4), either by changing the asthenospheric reference viscosity at the lithosphere base or by inserting a warm thermal anomaly simulating an ascending plume head (Fig. 2). We also test the influence of the lithosphere ductile strength that should modulate plate resistance to bending (2) by varying the mantle activation energy, $E_{\mathrm{a}}^{\mathrm{m}}$. At last, we further explore the TF mechanical structure (3) by imposing an increased width of the TF weak gouge, and different thermal structures of the plate boundary forming the TF.

2.5 Ranges of investigated physical properties

2.5.1 Brittle properties for oceanic crust, transform fault and mantle lithologies

The brittle parameter γ in Eq. (2) is related to the tectonic deviatoric stress, Δσ_xx, and to the lithostatic pressure, σ_zz (Turcotte and Schubert, 1982): Δσ_xx=γσ_zz. One may derive the relationship under compression between γ and the classical coefficient of static friction, f_s, defined by $f_{\mathrm{s}}=\mathit{\tau }/{\mathit{\sigma }}_{\mathrm{n}}$, where τ is the shear stress along the fault (Turcotte and Schubert, 1982):

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathit{\gamma }={\displaystyle \frac{\mathrm{2}{f}_{\mathrm{s}}(\mathrm{1}-\mathit{\lambda })(\mathrm{1}-{\mathit{\rho }}_{\mathrm{w}}/\mathit{\rho })}{\sqrt{\mathrm{1}+f_{\mathrm{s}}^{\mathrm{2}}}-f_{\mathrm{s}}}}\text{ if }{p}_{\mathrm{w}}\ne \mathrm{0}\mathrm{Pa},\text{(7)}& {\displaystyle }\mathit{\gamma }={\displaystyle \frac{\mathrm{2}{f}_{\mathrm{s}}}{\sqrt{\mathrm{1}+f_{\mathrm{s}}^{\mathrm{2}}}-f_{\mathrm{s}}}}\text{ if }{p}_{\mathrm{w}}=\mathrm{0}\mathrm{Pa},\end{array}$$

where λ is the pore fluid pressure coefficient, ρ_w is the water density, and p_w is the pore fluid pressure, assuming that p_w=ρ_wgz if λ=0 and p_w=ρgz if λ=1. The brittle parameter γ moderately depends on the average density in the overlying column, ρ (Fig. S1 in the Supplement). The internal friction coefficient, f_s, initially considered approximately constant (f_s∼0.6 to 0.85; Byerlee, 1978) is suggested to vary with composition from recent experimental data. For a dry basalt, f_s would be encompassed between 0.42 and 0.6 (Rocchi et al., 2003; Violay et al., 2012). Assuming high pore fluid pressure in the oceanic crust (λ≥0.45), γ_c from Eq. (6) is then close to 0.8 (Fig. S1). If the oceanic crust is altered by the formation of fibrous serpentine or lizardite, f_s decreases to 0.30 (Tesei et al., 2018), entailing γ_c∼0.05 if the pore fluid pressure is high (λ=0.9), which we consider the minimum realistic value for modeling the crustal brittle parameter (Fig. 3a). In the presence of chrysotile, f_s may even be reduced to 0.12 at low temperature and pressure (Moore et al., 2004), which would reduce γ_c to ∼0.01 (for λ=0.9), deemed as the extreme minimum value for γ_c. Note that relationship between the presence of fluid and its effect on the effective brittle strength (λ value) depends on the fault network and on the degree of pore connectivity, which may be highly variable (e.g., Carlson and Herrick, 1990; Tompkins and Christensen, 1999).

Regarding the TF, the fault material is assumed to be either mechanically similar to a weak serpentinized crust (γ_TF=0.05) or even softer (e.g., Behn et al., 2002; Hall and Gurnis, 2005). In that case, we set ${\mathit{\gamma }}_{\mathrm{TF}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$.

At mantle depths, the effect of pore fluid pressure on brittle strength is more questionable than at crustal levels. To simplify, we suppose the pore fluid pressure p_w to be very low, close to zero, assuming that the lithospheric mantle is dry in absence of any previous significant deformation. The coefficient of internal friction from Eq. (7) for a dry mantle decreases from f_s=0.65 (Byerlee, 1978) to f_s∼0.35 or 0.45 if peridotite is partly serpentinized (Raleigh and Paterson, 1965; Escartín et al., 1997), leading to γ_m between 2.8 and 0.8. However, assuming γ_m=2.8 would lead to an extremely high lithospheric strength (∼1 GPa at only 11 km depth) since our rheological model neglects other deformation mechanisms. We thus restrict the maximum γ_m to 1.6, which has been shown to allow for a realistic simulation of subduction force balance for steady-state subduction zones (Arcay et al., 2008). The most likely interval for γ_m is eventually [0.8–1.6] (Fig. 3b). The mantle brittle parameter γ_m might decrease to ∼0.15 (f_s=0.12) if chrysotile is stable, which is nevertheless unexpected at mantle conditions. Lower γ_m are considered unrealistic, even if γ_m=0.02 has been inferred to explain plate tectonic convection (in the case of a mantle devoid of a weak crustal layer; Korenaga, 2010).

Figure 3 Physical properties tested in this study and investigated ranges. (a) Brittle parameter for the oceanic crust, γ_c; (b) Brittle parameter for the mantle, γ_m; (c) Oceanic crust density, ρ_c; (d) Density of the weak medium forming the TF, ρ_TF; (e) Activation energy of the oceanic crust, $E_{\mathrm{a}}^{\mathrm{c}}$, assuming a non-Newtonian exponent n=3 in Eq. (4); (f) Lateral extent of the weak domain on both flanks of the TF, L_w. The parameter intervals vary from realistic ranges (in green) to extreme values (in yellow). They are still extended beyond these values, up to unrealistic ranges to achieve the conditions allowing for spontaneous subduction (in red).

2.5.2 Crust and transform fault densities

The oceanic crust density is varied from the classical value for a wet gabbro composition in the pressure–temperature conditions prevailing at the surface (2920 kg m⁻³; Bousquet et al., 1997; Tetreault and Buiter, 2014). Crust density in the blueschist facies reaches 3160 kg m⁻³, but we try even higher densities by imposing a mantle value. This would correspond to crust eclogitization and the heaviest crust to maximize the column weight within the older plate (OP) to promote its gravitational instability (Fig. 3c). Rocks forming the fault “gouge” are likely to be vertically highly variable in composition, possibly rich in buoyant phases such as serpentine and talc close to the surface (e.g., Cannat et al., 1991) and more depleted in hydrous phases at the deeper level. Below the Moho, down to its deepest portion, the fault may be compounded of a mix between oceanic crust and altered mantle (Cannat et al., 1991; Escartín and Cannat, 1999). The density of the fault gouge is thus likely to increase from the surface toward the deeper part of the fault, from a hydrated gabbro density to a mantle density. We thus test for ρ_TF values spanning from a gabbroic density to a mantle one (Fig. 3d). Note that these densities correspond to reference values at surface conditions (T=0 ^∘C and P=0 kbar), knowing that density is here a function of temperature through the coefficient of thermal expansion (Table 2).

2.5.3 Activation energy for the crust

The most realistic interval for the crustal activation energy $E_{\mathrm{a}}^{\mathrm{c}}$ can be defined from experimental estimates $E_{\mathrm{a}}^{\exp }$ for an oceanic crust composition. Nonetheless, $E_{\mathrm{a}}^{\exp }$ values are associated with specific power law exponent, n, in Eq. (4), while we prefer to keep n=3 in our numerical simulations for the sake of simplicity. Therefore, to infer the $E_{\mathrm{a}}^{\mathrm{c}}$ interval in our modeling using a non-Newtonian rheology, we assume that without external forcing, mantle flows will be comparable to sublithospheric mantle convective flows. The lithosphere thermal equilibrium obtained using a non-Newtonian rheology is equivalent to the one obtained with a Newtonian ductile law if the Newtonian E_a is equal to the non-Newtonian E_a multiplied by $\mathrm{2}/(n+\mathrm{1})$ (Dumoulin et al., 1999). As sublithospheric small-scale convection yields strain rates by the same order of plate tectonics ($\sim {\mathrm{10}}^{-\mathrm{14}}$ s⁻¹; Dumoulin et al., 1999), this relationship is used to rescale the activation energies experimentally measured in our numerical setup devoid of any external forcing. We hence compute the equivalent activation energy as follows: $E_{\mathrm{a}}^{\mathrm{c}}=(n+\mathrm{1})\times E_{\mathrm{a}}^{\exp }/(n_{\mathrm{e}}+\mathrm{1})$, where n_e is the experimentally defined power law exponent. The activation energy $E_{\mathrm{a}}^{\exp }$ in the dislocation creep regime is encompassed between the one for a microgabbro, 497 kJ mol⁻¹ (Wilks and Carter, 1990, with a non-Newtonian exponent n_e=3.4) and the one of a dry diabase, i.e., 485±30 kJ mol⁻¹ (Mackwell et al., 1998, with $n_{\mathrm{e}}=\mathrm{4.7}\pm \mathrm{0.6}$). For a basalt, $E_{\mathrm{a}}^{\exp }$ has been recently estimated to 456 kJ mol⁻¹ (Violay et al., 2012, with n_e∼3.6). Lower values inferred for other lithologies are possible but less likely, such as for a wet diorite ($E_{\mathrm{a}}^{\exp }=\mathrm{212}$ kJ mol⁻¹, n_e=2.4; Ranalli, 1995), and are used to define the lower bound of the “yellow” range for $E_{\mathrm{a}}^{\mathrm{c}}$ (Fig. 3e). A few experiments have shown that $E_{\mathrm{a}}^{\exp }$ can be as low as 132 kJ mol⁻¹ (n_e=3) if hornblende and plagioclase are present in high proportions (Yongsheng et al., 2009). This activation energy, as well as the one of a wet quartzite ($E_{\mathrm{a}}^{\exp }=\mathrm{154}$ kJ mol⁻¹, n_e=2.3; Ranalli, 1995), though used in numerous thermomechanical modeling studies of subduction, is considered an unrealistic value in a TF setting. Nevertheless, a low plate ductile strength promoted by a thick crust has been suggested to favor spontaneous subduction initiation at a passive margin (Nikolaeva et al., 2010). We choose to not vary the crustal thickness but to test in a set of experiments the effect of a very low crustal activation energy instead (equal to 185 kJ mol⁻³, Fig. 3e).

2.5.4 Distance from the transform fault of crust weakening

Regarding the lateral extent of the weak material, L_w, we test values in agreement with the observed large or relatively small TFs (L_w≤20 km, as described in the previous section) and increase them up to the extreme value of 50 km (Fig. 3f). The simulation results prompt us to perform experiments in which both lithospheres are entirely covered by the weak layer (L_w∼1110 km) to achieve the conditions of spontaneous subduction initiation.

2.6 Numerical code and resolution

The models are performed using the thermo-chemo-mechanical code of convection developed by Christensen (1992), which is based on an Eulerian and spline finite element method. Conservation equations are solved to obtain two scalar fields, which are temperature and stream function (Christensen, 1984). The simulation box is discretized into 407×119 nodes. The resolution is refined in x and z directions in the area encompassing the TF, i.e., between 966 and 1380 km away from the left-hand box side, and for depths shallower than 124 km, where node spacings are set to 1.67 km. Outside the refined domain, node spacing is 10.5 km in both directions. The tracer density is uniform over the simulation box (∼3.2 per km²), verifying that at least nine tracers fill the smallest meshes. This numerical discretization has been tested and validated in a previous study (Arcay, 2017). Note that because the total pressure is not directly solved by the code in Christensen (1992), the lithostatic pressure is used instead in Eq. (4).

The original code has been adapted to allow for the simulation of three different lithologies within the simulation box (Doin and Henry, 2001; Morency and Doin, 2004): the mantle, oceanic crust, and a weak layer that would mimic an altered or hydrated and, hence, weakened region around a TF, with specific densities and rheologies (see Sect. 2.3). Composition is tracked by markers advected along flow lines using a fourth-order Runge–Kutta scheme (van Keken et al., 1997).

3 Results

Here, we summarize first the experiments without OPS and then the simulations showing spontaneous gravitational instability of the OP. Next, we detail the effect of the different mechanical and geometrical parameters. Table 3 compiles the experiments explicitly quoted in the main paper. The exhaustive list of simulations performed in this study can be found in the Supplement: the experiments are compiled as a function of the plate age pair imposed at the TF in Table S1 in the Supplement, while they are ranked according to the simulated deformation regime in Table S2.

3.1 Overview of simulated behaviors other than old plate sinking

We obtain numerous behaviors different from OPS, varying as a function of (1) the plate age pair (A_y, A_o) and (2) the combination of densities, rheological parameters and the weak layer lateral extent (L_w). This large simulation set shown in Fig. 4 represents ∼73 % of the 302 experiments presented in this study, which do not show a clear OPS.

First, no tectonic deformation is modeled in many experiments, i.e., deformation only occurs within the asthenosphere below the plates but is almost totally absent at shallower depths where plate cooling takes place (Fig. 4.1). This is notably obtained if the YP is too old, that is, for A_y≥3 up to 17 Myr depending on the physical parameter set (Fig. 6).

Second, we observe the YP ductile dripping, leading to the plate dismantlement, corresponding to a series of several fast lithospheric drips, soon after the simulation start (Fig. 4.2), modeled when ductile strengths are low. The OP is not affected and solely cools through time.

Third, a transient retreat of the YP is modeled, in very few experiments, while the OP remains motionless (Fig. 4.3). This occurs if the YP is very young (A_y≤2 Myr) and if the TF density, ρ_TF, is low (equal to the gabbro density). Because of its buoyancy, the weak material forming the TF rises up to the surface as soon as simulation starts. This fast vertical motion (velocities ≥50 cm yr⁻¹) is partly transmitted horizontally and deforms the weaker and younger plate, triggering a backward motion. Velocities vanish as plate cooling proceeds.

Fourth, the YP sinking is triggered in some models (Fig. 4.4). The gravitational instability of the YP is very similar to the one expected for a thick plate spontaneous sinking (as sketched in Fig. 1). The polarity of the YP sinking depends on the density imposed for the TF interface (ρ_TF) and on whether (or not) the weak layer covers the YP surface (L_w>0 km). The duration of YP spontaneous sinking is always very brief (<0.5 Myr): either the process does not propagate fast enough to compete against plate cooling and strengthening (Fig. 4.4b) or the diving YP segment is limited by the imposed length L_w of lithosphere recovered by the weak material (Fig. 4.4a).

Fifth, in one experiment, a double subduction initiation is observed: while YP sinking initiates, the OP also becomes unstable and starts sinking when a wider portion of weak and dense material (L_w=50 km) is included (Fig. 4.5). Nevertheless, the OP slab rapidly undergoes slab break off once the L_w long weak segment has been entirely subducted (Fig. 4.5, 0.62 Myr), which we deem as too short to represent a successful OPS initiation since the subducted slab length is limited to 50 km.

Sixth, the vertical subduction of the YP initiates at the TF when the TF material is as dense as the mantle and vertically drags the YP into the mantle (Fig. 4.6). The motion can be transmitted away from the TF up to 500 km backward but systematically entails a YP stretching at the surface, as the slab is young and soft (A_y≤7 Myr). This prevents subduction from lasting more than 1.5 Myr. Moreover, plate cooling frequently freezes the downward YP flow (Fig. 5.6, bottom row).

Finally, in ∼40 % of experiments in which OPS initiation appears to start, the process freezes up and does not evolve into a developed plate sinking. The OP bending stops very early, typically in less than ∼0.4 Myr (Fig. 4.7), especially for OP older than 80 Myr. The velocities within the OP then vanish quite fast (Fig. 4.7a). OPS also aborts even when the mechanical decoupling does occur at the TF if hot mantle flows are too slow and/or if the lateral extent of the weak material L_w is narrow (Fig. 4.7b).

3.2 Modes of OPS triggering

Spontaneous subduction is modeled when one of the two lithospheres is gravitationally unstable, which occurs if the total lateral density offset (vertically integrated) at the plate boundary is not balanced by plate, mantle, and TF resistance to deformation, as summarized in Sect. 2.4. We observe the spontaneous sinking of the OP for quite various pairs of lithosphere ages (Fig. 5), which mostly depends on the chosen set of rheological parameters and on the presence of the weak layer at the whole plate surface. When simulated, OPS occurs following one of two basic ways, later called mode 1 and mode 2. Mode 1 happens in approximately one-half of OPS cases (Fig. 5a) and is the closest to the mechanism envisioned in the spontaneous subduction concept (Fig. 1b). The mantle flow generated by the OP sinking triggers an asthenospheric upwelling focusing along the weak TF “channel” up to the surface (“asthenosphere invasion” in Fig. 1b), while the YP remains mostly motionless. The subduction process develops due to a fast hinge rollback. As mantle velocities are huge, exceeding tens of meters per year in many cases, the asthenosphere catastrophically invades the box surface, filling a domain that is soon larger than 200 km, as depicted in Fig. 5a.

In mode 2, asthenosphere invasion does not occur at the surface and is often limited to the YP Moho. Mantle flow induced by OP bending drags the YP toward the OP (Fig. 5b, c). As a consequence, a significant mass of dense crust is transferred from the top of the YP to the one of the OP, where the accumulated crust builds a crustal prism that loads the OP, amplifying its bending and sinking. This phenomenon is observed in numerous cases, systematically if the YP age is 2 Myr (Table 3), and in several cases when A_y is either 0 or 5 Myr (simulations S1a to S2b, S22j–k). In both initiation modes, velocities at the slab extremity are very high (14.6 cm yr⁻¹ in simulation S1a, 0 vs. 2, up to ∼ 180 cm yr⁻¹ in simulations S10a, 0 vs. 80 and S11a, 0 vs. 100). The duration to form a slab longer than ∼200 km is less than 1.5 Myr. The kinetics of the OPS process modeled in this study are consequently always very fast. This swiftness most likely comes from the significant weakness that must be imposed in our modeling setup to obtain OPS triggering (see Sect. 3.3.2).

Figure 4 Illustration of the different simulated behaviors, OPS apart: close-up on the transform fault. (1) Absence of plate deformation (simulation S37x, Table 3). (2) Young plate dripping and dismantlement (simulation S17f). (3) YP retreat (simulation S16c). (4) Initiation of YP transient sinking (simulation S16b, panel a, and simulation S36b, panel b). (5) Simultaneous initiation of YP and OP sinking processes (simulation S14n). (6) Initiation of YP vertical subduction at the TF (simulation S17o). (7) OP sinking initiation that soon aborts (simulation S33a, panel a, and simulation S16a, panel b). No vertical exaggeration. The velocity scale depicted in green is specific to each simulation. The parameter boxes are color coded as a function of the investigated ranges depicted in Fig. 3.

Figure 5 Illustration of OPS: mode 1 in simulation S27c (panel a) vs. mode 2 in simulations S14i (panel b), and S22j (panel c). No vertical exaggeration. The parameter boxes are color coded as a function of the investigated ranges depicted in Fig. 3. Note that the velocity scale in panel (c) is specific for each snapshot. The dashed lines in the middle sketch are a schematic outline of the stream function, and the green arrows illustrate velocities.

3.3 Influence of tested parameters

The regime diagrams displayed as a function of the plate age pair (A_y,A_o) sum up our main results obtained as a function of the assumed rheological set, density field, and the lithological distribution at the surface (oceanic crust vs. TF weak material; Fig. 6). These eight regime diagrams bring out the respective influence of the main physical parameters tested in this paper, especially for deciphering conditions allowing for OPS. YP dismantlement, basically occurring when the ductile crust is softened, is not represented in the regime diagrams (discussed at the end of the section).

Figure 6 Regime diagrams as a function of the combination of rheological properties, densities, and TF weak domain extent, L_w. “Irrelevant domain” would correspond to cases where A_y>A_o. In dotted areas, only transient behaviors, lasting less than 0.5 Myr, are modeled. Inside the inserted sketches, the numbers refer to the panel numbering in Fig. 4. The box colors correspond to parameter ranges depicted in Fig. 3. The parameter combination γ_c, γ_m, ρ_c, ρ_TF, and L_w, is, respectively, in panel (a): 0.05, 1.6, 2920, 2920 kg m⁻³, and 0 km; in panel (b) 0.05, 1.6, 2920 kg m⁻³, 3300, and 0 km; in panel (c): 0.05, 1.6, 2920, 3300 kg m⁻³, and 50 km; in panel (d): 0.05, 1.6, 3300, 3300 kg m⁻³, and 0 km; in panel (e): 0.05, 1.6, 3300, 3300 kg m⁻³, and 50 km; in panel (f): $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$, 1.6, 3300, 3300 kg m⁻³, and 1100 km; in panel (g): $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$, 0.1, 3300, 3300 kg m⁻³, and 1100 km; and in panel (h): $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$, 0.05, 3300, 3300 kg m⁻³, and 1100 km. In all panels, $E_{\mathrm{a}}^{\mathrm{c}}=\mathrm{360}$ kJ mol⁻¹. The corresponding experiments are displayed in Fig. S2 and S3 in the Supplement. In panel (f), the boundary between the “no subduction” and “OPS” domains corresponds to the relationship $A_{\mathrm{o}}/A_{\mathrm{y}}^{\mathrm{2.5}}\mathit{\gtrsim }\mathrm{0.75}$ Myr^−1.5. When OPS is simulated (panels e to h), the conditions in A_o-A_y prevailing at subduction initiation inferred for Yap, IBM, and Matthew and Hunter (Table 1) are superimposed on the regime diagrams.

3.3.1 Transform fault and oceanic crust densities

Densities strongly affect the evolution of the TF system. If the TF weak medium is buoyant (${\mathit{\rho }}_{\mathrm{TF}}={\mathit{\rho }}_{\mathrm{c}}=\mathrm{2920}$ kg m⁻³), the TF material rises up to the surface forming a small and localized buoyant diapir that pushes laterally on the younger lithosphere (Fig. 4.3). The YP either shortens if it is weak enough (A_y≤2 Myr, Fig. 6a) in a backward motion or starts sinking if the YP thickness is intermediate ($\mathrm{2}<A_{\mathrm{y}}<\mathrm{20}$ Myr). On the other hand, a heavy material filling the TF gouge (ρ_TF=3300 kg m⁻³) inverts the aforementioned mechanics by pulling the YP downward at the TF to form a vertical subduction (Fig. 4.6, labeled YPVSI for “YP vertical subduction initiation” in Table 3). Note that when the fault density ρ_TF is very high, the oceanic crust density, ρ_c, buoyant or not, does not actually affect the mode of YP deformation (compare diagrams b and d in Fig. 6).

3.3.2 Lateral extent of the weak material

The results presented in Sect. 3.3.1 are obtained when the weak material is localized at the TF only (L_w=0 km). Assuming that the weak material laterally spreads out away from the TF (L_w>0 km), the mode of YP vertical subduction switches to YP sinking by gravitational instability. This is observed when young plates are modeled on both sides of the TF (A_y<5 Myr, A_o<40 Myr, Fig. 6c). The boundary between the dense weak material and the buoyant and stronger oceanic crust more or less acts as a “secondary” plate boundary, decoupling the two lithological parts of the YP, which does not occur if there is no buoyancy contrast between the crust and the weak material (Fig. 6e).

Moreover, we observe that enlarging the weak domain enables OPS in some cases if the YP is very thin (A_y≤2 Myr), regardless of the oceanic crust density (Fig. 6c, e), although OPS aborts fast, as OP subduction is limited to the weakened length L_w (set to 50 km, Fig. 6e). Simulations show that OP sinking is enhanced if L_w is much wider than expected in nature (L_w≥50 km, Fig. 3f, g, h). Otherwise, the backward propagation of bending is hindered, which stops the OPS process. We conclude that a very wide area of crust weakening on both sides of the TF is a necessary condition to simulate OPS. We quantify more accurately for different pairs of plate ages with minimum length L_w, allowing for a developed OPS in the Supplement (Sect. S2). These age pairs are selected to cover a wide range of YP ages (2 to 20 Myr). We find that the domain of weakened crust to impose in the vicinity of the TF is too large to be realistic, at least for classical mantle rheology, with the only exception being the setting with a very thin YP (A_y=2 Myr). These results suggest the strong resistant characteristic of thick YP in OPS triggering.

3.3.3 Crust brittle strength

What is the threshold in crust weakening enabling OPS? A usual value of the crust brittle parameter (γ_c=0.05) does not allow for OPS (Fig. 6a to e). Our simulations show that if γ_c is 100 times lower (${\mathit{\gamma }}_{\mathrm{c}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$), OPS can initiate for numerous plate age pairs if the whole crust is mechanically weak (L_w=1100 km, Fig. 6f) but such a brittle parameter seems unrealistic. To determine the threshold in γ_c allowing for OPS, we choose a high plate age offset, 2 vs. 80, the most propitious for OPS (keeping L_w=1100 km). We determine that the threshold in γ_c is encompassed between 10⁻³ and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ (simulations S18b, c, d, and e), which is still less than the lower bound of acceptable γ_c ranges (Fig. 3a). We hypothesize that for a small plate age offset, the threshold in γ_c would have to be even lower to observe OPS triggering.

3.3.4 Plate bending and mantle brittle parameter

Surprisingly, a very low crust brittle parameter is not sufficient for simulating OPS for some large plate age offsets, such as for (A_y,A_o) = 10 vs. 100 or 5 vs. 120 (simulations S41a and S29b, Table 3, Fig. 6f). A mechanism is thus hindering OPS. We assume that thick OPs are too strong to allow for bending. We test this by reducing the mantle brittle parameter, γ_m, that affects the maximum lithospheric stress in our brittle–viscous rheology, from 1.6 (Fig. 6f) to 0.1 (Fig. 6g) and 0.05 (Fig. 6h). The domain of the plate age pair where OPS can occur is then greatly enlarged toward much lower plate age offsets. We note that in most experiments showing “mantle weakening”-induced OPS, OPS stops by an early slab break off, once the infant slab reaches 200 to 300 km length because the reduced slab strength cannot sustain a significant slab pull (Fig. 5c).

In a limited set of experiments, we determine the threshold in ${\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{crit}}$ below which OPS occurs. This threshold depends on the OP age: ${\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{crit}}\sim \mathrm{0.06}$ for the plate age pair 10 vs. 40 (simulations S37c to f) but the threshold is higher (${\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{crit}}\ge \mathrm{0.1}$) for the plate age pair 10 vs. 50 (sim. S38a). The thicker the OP is, the easier the OPS triggering, as one may expect. We next compare experiments in which the OP and the YPs are both progressively thickened by considering the following age pairs: 10 vs. 50 (sim. S38a), 15 vs. 60 (S47a), 20 vs. 80 (S51a–e), and 25 vs. 100 (S57a–b). The experiments show that ${\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{crit}}$ is ≥0.1, ≥0.1, ∼0.07, and ∼0.06, respectively. Hence, the plate rigidity has to be reduced as YP thickness increases, despite the joint OP thickening, down to extremely weak γ_m ranges (${\mathit{\gamma }}_{\mathrm{m}}^{\mathrm{crit}}\ll \mathrm{0.1}$, Fig. 3). Despite the driving influence of thicker OP, thickening the YP impedes OPS in a much stronger way. Moreover, we test different means to lower the OP rigidity. For four plate age pairs for which OPS aborts (5 vs. 35, 7 vs. 70, 7 vs. 80, and 7 vs. 90), we decrease the mantle ductile strength by lowering the activation energy $E_{\mathrm{a}}^{\mathrm{m}}$ (Table 2) but keep constant the mantle viscosity at 100 km depth and the mantle brittle parameter (γ_m=1.6). We find that lowering $E_{\mathrm{a}}^{\mathrm{m}}$ instead of the mantle brittle parameter is much more inefficient for obtaining OPS (Table S1). Finally, our results suggest that the factors of the resistance to OPS mainly come from OP flexural rigidity and YP thickness and stiffness, in agreement with previous studies (Nikolaeva et al., 2010; Zhou et al., 2018).

3.3.5 Ductile strength decrease

The main effect of imposing a decrease in the crust and TF ductile strength (lowering $E_{\mathrm{a}}^{\mathrm{c}}$ to 185 kJ mol⁻¹) is to trigger the fast dismantlement of YP by lithosphere dripping if the YP is young (A_y=2 Myr, Fig. 4.2). Otherwise, a low $E_{\mathrm{a}}^{\mathrm{c}}$ has no effect on the deformation of the two plates. One exception appears in simulation S14e, in which the weak ductile strength triggers mode 2 OPS. In this particular setup, both lithospheres are very thin ($(A_{\mathrm{y}},A_{\mathrm{o}})=\mathrm{2}$ vs. 5) and could be regarded as “crustal” plates because the mantle lithosphere is very thin or almost absent. In this simulation, the YP strength profile is actually similar to the other cases yielding mode 2 OPS (see Sect. S4 in the Supplement), which should explain why decreasing $E_{\mathrm{a}}^{\mathrm{c}}$ allows for OPS in this unique case. The YP destabilization and dripping result from the high crust density (ρ_c=3300 kg m⁻³) assumed in experiments performed with a reduced $E_{\mathrm{a}}^{\mathrm{c}}$. Indeed, in experiments using a usual crustal density (ρ_c=2920 kg m⁻³), YP vertical subduction is obtained instead (simulation S15g, 2 vs. 10, for instance).

3.3.6 Plume-like thermal anomaly

The thermal anomaly simulating an ascending plume head below the TF produces effects very similar to those of a reduced $E_{\mathrm{a}}^{\mathrm{c}}$: no effect if plates are older than 2 Myr and YP dismantlement if A_y=2 Myr and if the crust is dense (ρ_c=3300 kg m⁻³). Otherwise, for a normal crust density, a short stage of YP vertical subduction occurs after plume impact (2 vs. 10, simulation S15h). The hot thermal anomaly never triggers OPS in our modeling, contrary to other studies, even if we have investigated large plate age contrasts (2 vs. 40, sim. S17j, and 2 vs. 80, S18k) as well as small age offsets and plates younger than 15 Myr (Table S1). To obtain a successful plume-induced subduction initiation, it has been shown that the plume buoyancy has to exceed the local lithospheric (plastic) strength. This condition is reached either when the lithosphere friction coefficient is lower than ∼0.1 (Crameri and Tackley, 2016), and/or when the impacted lithosphere is younger than 15 Myr (Ueda et al., 2008), or when a significant magmatism-related weakening is implemented (Ueda et al., 2008) or assumed (Baes et al., 2016) in experiments reproducing modern Earth conditions. We hypothesize that if the mantle brittle parameter was sufficiently decreased, we would also achieve OPS by plume head impact. Besides, lithosphere fragmentation is observed by Ueda et al. (2008) when the plume size is relatively large in relation to the lithosphere thickness, in agreement with our simulation results showing the dismantlement for a significantly young (A_y=2 Myr) and thin lithosphere.

3.3.7 Additional tests on OPS conditions: fault gouge strength and width, transform fault vs. fracture zones, and asthenosphere viscosity

We sum up in this section the extra experiments performed to provide the precision of the mechanisms involved in OPS triggering. The detailed results are described in Sect. S3 in the Supplement. We first test the necessity of the fault softness to simulate OPS by inverting the oceanic crust and TF respective brittle parameter for models that originally displayed OPS (thus by setting for the inversion experiments: γ_TF=0.05, while γ_c=0.0005). We find that a very low TF strength is critical to model OPS.

We next wonder if OPS (when not modeled) could be triggered by widening the fault gouge from the surface to the bottom of the fault (domain 1 in Fig. 2) by setting the fault width to 20 km instead of 8.3 km in experiments that did not initially show OPS. The simulations show that OPS still does not occur, even if the mechanical decoupling is maximized (γ_TF decreased to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$). The mechanical interplate decoupling is hence not sufficient alone to trigger OPS, at least for a 20 km wide fault.

Subsequently, we investigate the possible role of the TF thermal structure. The interplate domain is assumed to be very thin and modeled by a stair-step function (Sect. 2.2). In nature, this setup would correspond to an active transform fault. If the fault is instead inactive (fracture zone), the thermal state of the plate boundary is likely to be cooled by thermal conduction, possibly stronger and more resistant to plate decoupling. We test the effect of the TF thermal structure for the two plate age pairs for which OPS is simulated when crust weakening is assumed (γ_c=0.0005) by widening the TF thermal transition (from 11 up to 70 km), keeping the weak material forming the fault gouge at the center of the thermal transition in all cases. All these experiments show OPS. We here verify that the fault gouge weakening, governed by the soft material brittle properties, is independent of temperature and, at first order, is independent of the fault activity in our 2-D setup.

We finally test if a decrease in asthenosphere strength could help OPS triggering by unbalancing the OP weight excess. The experiments show that asthenospheric velocities and OP deformation are slightly amplified but still not enough to trigger OPS.

Table 4 Relationship between investigated thermomechanical parameters and the main forces involved in the spontaneous subduction initiation modeled in this study. When increased from an experiment to another, the mechanical parameters are interpreted as either favoring OPS (plus sign), hampering OPS (minus sign), or not affecting OPS triggering (zero). Empty cells indicate that the relation between parameter and force cannot be inferred from our results, while a question mark means that the assumed relationship is not clear.

4 Analysis: mechanism and conditions of OPS triggering

We derive the main processes involved in successful OPS triggering from the results presented in Sect. 3. In the following paragraphs, we discuss the relationship between the different thermomechanical parameters and their influence on the forces driving and resisting to OPS, which are summed up in Table 4.

4.1 OPS mode 1 vs. mode 2

The velocity fields show that modes 1 and 2 differ by the resulting YP kinematics. When mode 1 occurs, the YP remains almost motionless with respect to underlying asthenospheric flows (speeds ≤1 and ≥20 cm yr⁻¹, respectively; see Supplement Fig. S7). In contrast, velocities during mode 2 are closer between the YP and the asthenosphere, where speeds are high (between 25 and 100 cm yr⁻¹). Moreover, within the set of simulations showing OPS (Sect. 3.2), we find that mode 2 occurs in simulations where the YP age is 2 Myr for various rheological sets, or if A_y=5 Myr, provided that the mantle brittle strength is reduced. In all these experiments, the strength of the YP bottom part is the closest to the asthenospheric one (viscosity ratios ≲10² to 10³; Fig. S8 in the Supplement). In contrast, the focusing of asthenosphere flows toward the weak TF is observed when the viscosity offset between the YP and the underlying mantle exceeds 10². We hypothesize that mode 2 results from a strong coupling between the YP and the asthenosphere, which is related to the asthenosphere ductile strength (Table 4). The particular (though not meaningful) case of A_y=0 Myr is addressed in the Supplement (Sect. S4).

The mode 2 OPS may be envisioned as an asymmetric double-sided subduction (Gerya et al., 2008, see the sketch in Fig. 5). In this subduction mode, the thick OP sinking drives the YP downward flow at the proto-slab surface because plate decoupling at shallow levels does not occur. The shallow interplate decoupling is hence not required in mode 2 since the YP is easily dragged by asthenosphere. In contrast, asthenospheric flows related to OPS in mode 1 might not be able to drag the YP because of the high viscosity offset between YP and asthenosphere . The asthenosphere upwelling along the TF (“upwelling force” in Table 4) would then result from the need to decouple the respective motions of the two plates, i.e., to accommodate OP downwelling and hinge retreat, whereas the YP is almost motionless. Simulations by Gerya et al. (2008) suggest that the thick OP lubrication by metasomatism is a way to force plate decoupling to model one-sided subduction.

4.2 Plate ages allowing for old plate sinking

The boundary between OPS and the absence of subduction can be defined for a normal mantle brittle strength γ_m=1.6 (Fig. 6f) using simulations in which OPS aborts (such as simulations S25a, 5 vs. 35; S29b, 5 vs. 120, or S33a, 7 vs. 80, Fig. S3 in the Supplement). We observe a dichotomy in the OPS domain boundaries. On the one hand, for thick OP (A_o>100 Myr), OPS is prevented if the YP is not extremely thin (plate age younger than 5 Myr). On the other hand, for thinner OP (A_o≤100 Myr), we experimentally show that the OPS condition corresponds to the following relationship: $A_{\mathrm{o}}/A_{\mathrm{y}}^{\mathrm{2.5}}\mathit{\gtrsim }\mathrm{0.75}$ Myr^−1.5 (Fig. 6f). In both cases, the influence of A_y is either strong or predominant. The YP age is the major determining factor in the TF evolution, compared to the OP age (separately considering the cases where the mantle brittle strength is reduced), which confirms the conclusion derived in Sect. 3.3.4 on the highly resistant effect of the YP thickness. This hindering effect results from two processes. On the one hand, high A_y ages yield low-pressure gradients across the TF due to a density contrast that decreases with YP aging (e.g., Hall and Gurnis, 2003). On the other hand, YP aging increases the YP strength competing against asthenosphere upwelling in the vicinity of the TF in OPS mode 1 and YP stretching far away from the TF to accommodate YP dragging in mode 2 (Table 4). As a result, the conditions that are the most propitious for OPS correspond to TFs, where the thinner lithosphere is as young as possible.

4.3 Parameters resisting and promoting old plate sinking

OPS is triggered if the pressure gradients at the TF related to the density offset exceed plate and mantle resistance to deformation. Density contrasts are maximized when the YP is thin, which partly explains the dominant role of the age of the YP (“proto-overriding” plate), compared to A_o, on subduction initiation, as already underlined in other studies (Nikolaeva et al., 2010; Dymkova and Gerya, 2013). Our results show that plate instability is essentially promoted by three mechanical conditions: when low brittle strengths are assumed for (1) the oceanic crust and (2) the mantle and (3) if the TF allows for plate decoupling. A weak brittle crust (1) enhances the fast propagation of deformation at shallow depths, which cannot be obtained in our modeling by crust ductile softening (contrary to Nikolaeva et al., 2010, in a passive margin setup). Moreover, the lowering of the crust brittle strength must be developed far away from the TF to allow for OPS. Although the minimum spatial scale of crust softening depends on plate age pair, we find that it is generally on the order of a hundred(s) of kilometers. Low brittle mantle strength (2) strongly promotes not only the OP plate bending and sinking by limiting the plate flexural rigidity but also YP deformations close to the plate boundary where asthenosphere upwelling focuses in mode 1. Finally, the TF must also be weak to enable mechanical decoupling between neighboring plates (3).

The amplitude of the preceding processes is regulated by five of the six physical parameters investigated in this study, as the activation energy $E_{\mathrm{a}}^{\mathrm{c}}$ does not actually affect OPS triggering (Fig. 3). Clearly, OPS cannot be simulated for a realistic set of physical parameters (Fig. 6a). To achieve OPS, the cursors controlling the plate mechanical structures have been tuned beyond the most realistic ranges (“yellow” domain, Fig. 3) for two parameters at least, and beyond reasonable values for at least one parameter (“red” domain, Fig. 6e to h). Nevertheless, combining different unlikely (“yellow”) parameter values (for ρ_TF and L_w) does help to achieve OPS for slightly less extreme mechanical conditions, as one parameter only has to be pushed up to the unrealistic (“red”) range (ρ_c, Fig. 6e). Note, however, that the plate age intervals showing OPS are then extremely narrow (A_y<3 Myr, A_o<25 Myr) and are not consistent with the three potential candidates of natural OPS.

5 Discussion

5.1 Model limitations

5.1.1 2-D versus 3-D setup

Subduction initiation at a TF is here simplified using a 2-D process, whereas the fault strike-slip basically implies that it takes place as a 3-D phenomenon. For instance, the setting of Matthew and Hunter used to be a subduction-transform edge propagator (“STEP”) fault 2 Myr ago (Patriat et al., 2015), where 3-D mantle flows associated with the Australian plate subduction likely affected the TF structure and evolution, possibly favoring subduction initiation. This case exemplifies the role of 3-D far-field tectonics during subduction infancy (Table 1) and the potential role of deep mantle flows. Upward and downward mantle flows, even far away from the initiation site, have been shown to be able to initiate subduction in 2-D models (Lu et al., 2015; Baes and Sobolev, 2017). On the other hand, studies in 3-D by Boutelier and Beckett (2018) and Zhou et al. (2018) showed that subduction initiation depends on along-strike variations in plate structure. However, the strike-slip kinematics of an active TF have up to now, to our knowledge, never been taken into account in subduction initiation simulations. We show that a TF thermal state cooler than that modeled by a stair-step function does not hinder OPS. However, our results verify that in 2-D, without simulating the TF strike-slip, the process of spontaneous OPS has to occur at simulation onset to avoid the impeding effect of plate stiffening with further conductive cooling. Including the TF motion in 3-D experiments would compete against this strengthening effect in the area nearby the active spreading center.

Finally, one may argue that a 3-D setup would intrinsically facilitate OPS propagation at a transform fault. Plate sinking might initiate at the location where the offset in plate thickness is maximum (in the vicinity of a ridge spreading center) and then propagate away from this point (Zhou et al., 2018). However, as we focus on subduction initiation strictly speaking and not on subduction propagation, the use of a 2-D setup should remain meaningful to unravel the conditions of spontaneous sinking for a given plate age pair, considering apart the problem of the transform fault slip.

5.1.2 Free slip vs. free surface condition

Our results show that in most of our experiments showing OPS, the oceanic crust must be significantly weak (Fig. 6f–g). One may argue that the necessity of decoupling propagation close to the surface by shallow softening is related in our modeling to the absence of free surface (e.g., Crameri and Tackley, 2016). We test it by seeking for the threshold in the crustal brittle parameter allowing for OPS for one plate age pair 5 vs. 40 (sim. S26a in Table 3) as a function of the mechanical boundary condition imposed at the box top, either free-slip without vertical motion or free surface, mimicked by inserting a “sticky water” layer (see the Supplement Sect. S3 and Fig. S6). For the selected plate age pair, the threshold in crustal brittle parameter turns out to increase from ∼0.0025 without free surface to ∼0.0175. Hence, the necessary crust weakness that must be imposed to model OPS may be overestimated by a factor ∼7. This result agrees with previous studies showing that the free surface condition promotes the triggering of one-sided subduction in global mantle convection models (Crameri et al., 2012). Nevertheless, note that the threshold enabling OPS when the free surface is taken into account is still an unlikely value, since it is close to the limit of the extremely low range of the crust brittle parameter (“red” domain, Fig. 3).

5.1.3 Initiation swiftness and influence of elastic rheology

In a TF or fracture zone numerical setting without any external forcing, if subduction initiation has to occur, it can only take place at simulation onset because plate cooling results in a fast stiffening of oceanic lithospheres and, second, quickly attenuates the plate density offset (Hall and Gurnis, 2003). The process of subduction initiation modeled in our study systematically occurs very briefly after the simulation start, in less than 1 to 1.5 Myr. This quite “catastrophic” way of initiation has also been simulated in less than 0.8 Myr for other tectonic settings or triggering modes, such as passive margins (Nikolaeva et al., 2010; Marques and Kaus, 2016) or plume-induced mantle flows (Lu et al., 2015), using rheological conditions very similar to the ones assumed in this study. The initiation process is slightly slowed down but remains fast (duration <3 Myr) when the necessary weakness of the plate's stronger part is not fully imposed at simulation onset but progressively develops due to damaging or water-related weakening effects (Hall and Gurnis, 2003; Gurnis et al., 2004; Gerya et al., 2008; Dymkova and Gerya, 2013). Moreover, such unrealistically high velocities at plate sinking onset may result at least in part from the 2-D setup since, in a 3-D setup, the along-strike propagation slows down the initiation process; however, speeds of hinge retreat remain significantly high (between 13 and 20 cm yr⁻¹ in Zhou et al., 2018). In addition, by neglecting elastic deformation, the amount of plate and interplate weakening required to trigger OPS may be excessive (Farrington et al., 2014). Nonetheless, the potential effect of elasticity on the OPS kinetics is not clear. On the one hand, including elasticity could slow down OPS initiation by increasing the threshold in the strength contrast, as aforementioned. On the other hand, the incipient subduction has been shown to remain as fast as modeled in the present study in elasto-visco-plastic models testing different modes of subduction initiation (Hall and Gurnis, 2003; Thielmann and Kaus, 2012; Baes et al., 2016). Previous modeling of subduction initiation including elasticity showed that the elastic flexure was a basic term of the subduction force balance (McKenzie, 1977; Hall and Gurnis, 2003; Gurnis et al., 2004). In our model, plate bending occurs by viscous (and brittle) deformations, as in numerous approaches to subduction that have successfully reproduced topographies and the strain and stress patterns observed in natural cases (e.g., Billen and Gurnis, 2005; Buffett, 2006). However, if elasticity might compete against subduction initiation by limiting the localization of lithospheric shearing, it may also help incipient subduction through the following release of stored elastic work (Thielmann and Kaus, 2012; Crameri and Tackley, 2016). Consequently, the threshold in mechanical parameters necessary to achieve OPS would probably be offset if elasticity was included.

5.1.4 Weakening of the oceanic mantle lithosphere

Our results show that OPS is strongly facilitated if the lithospheric mantle is weak. The reduction in mantle strength depends on the plate age pair and on other mechanical parameters, such as γ_m. A first-order estimate of the necessary mantle weakening is computed by comparing cases showing OPS to those in which OPS fails (Sect. S5 in the Supplement). The mantle weakening allowing for OPS is low to moderate for young plates and high plate age offsets (strength ratio ≤35), and larger when the plate age contrast is small (strength ratio ∼280). One may wonder if such mantle strength decreases are realistic. Different mechanisms of mantle weakening that we do not model may be discussed, such as (1) low-temperature plasticity (Goetze and Evans, 1979), that enhances the deformation of slab and plate base (Garel et al., 2014); (2) dislocation-accommodated grain-boundary sliding (GBS; Hansen et al., 2012); (3) grain-size reduction when diffusion linear creep is activated; or (4) fluid-related weakening. Peierls plasticity (1) limits the ductile strength in a high-stress regime at moderately high temperatures (≲1000 ^∘C; Demouchy et al., 2013) but requires a high differential stress (>100 to 200 MPa) to be activated. Similarly, GBS power law regime (2) operates if stresses are >100 MPa, for large strain and low temperature (<800 ^∘C; Drury, 2005). In our experiments, the simulated deviatoric stress is generally much lower than 100 MPa (Sect. S5 in the Supplement). For such low stresses, neither the Peierls nor the GBS creeps would be activated, hence we do not expect a major change in our results if they were implemented. Grain-size sensitive (GSS) diffusion linear creep (3) can strongly localize deformation at high temperature (e.g., Karato et al., 1986). In nature, GSS creep has been observed in mantle shear zones in the vicinity of a fossil ridge in Oman in contrast at rather low temperature (≲1000 ^∘C; Michibayashi and Mainprice, 2004), forming very narrow shear zones (<1 km wide). However, the observed grain-size reduction of olivine is limited to ∼0.2–0.7 mm, which cannot result in a noticeable viscosity reduction. A significant strength decrease associated with GSS linear creep requires additional fluid percolation once shear localization is well developed within the subcontinental mantle (e.g., Hidas et al., 2016). The origin of such fluids at great depth within an oceanic young lithosphere is not obvious. Furthermore, GSS linear creep may only operate at stresses <10 MPa (Burov, 2011), which is not verified in our simulations (Sect. S5 in the Supplement).

In our models, mantle weakening is achieved by decreasing the mantle brittle parameter γ_m, to mimic the weakening effect of hydrated phases (4), such as talc or serpentine minerals (Sect. 2.5.1). Dymkova and Gerya (2013) show that the percolation of sea water down to ∼25 km depth during early OP deformation can enable the thick plate bending, assuming low porosity (≤2.5 %) and low mantle matrix permeability (10⁻²¹ m²) to significantly increase pore fluid pressure. In our approach, a high pore fluid pressure ratio (λ>0.5) is associated with a low mantle brittle parameter (γ_m<1, Fig. S1), for which OPS is modeled for a broad ranges of plate ages (Fig. 6g–h), in agreement with results by Dymkova and Gerya (2013). However, the low permeabilities assumed by Dymkova and Gerya (2013) are questioned by recent experiments of mantle hydration at ridge and of water percolation in a peridotite, and by estimates from a peridotite aquifer (Dewandel et al., 2004; Godard et al., 2013; Farough et al., 2016). These studies rather infer permeabilities between 10⁻¹⁹ and 10⁻¹⁶ m², which would hamper high pore fluid pressures and, eventually, plate bending.

5.2 Comparison between OPS model requirements and natural cases

When analyzing the results of our parametric study (Sect. 4), the striking conclusion is that none of the realistic sets of parameters allowed for spontaneous subduction. Figure 6e, f, g, and h show that, even if one of the plates is extremely young (<7 Myr), the oceanic crust should be very dense (ρ_c=3300 kg m⁻³), as well as drastically weakened (${\mathit{\gamma }}_{\mathrm{c}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$) at considerable distances from the TF (L_w≥50 km), to satisfy OPS necessary conditions. Assuming that such rather extreme conditions were fulfilled, OPS must develop (1) at simulation onset before plate cooling so that gravitational instability is maximal and (2) catastrophically in terms of the kinetics of the process, with sinking rates ≥15 up to 180 cm yr⁻¹.

As depicted in Sect. 1, only two natural cases, IBM and Yap, attest to subduction initiation of an old oceanic plate beneath a young one at a TF, which later evolved into a mature subduction during the Cenozoic. The range of ages for both plates at time of initiation (Table 1) has been plotted in Fig. 6f, g and h. The plate ages at Yap subduction initiation are incompatible with the conditions of OPS inferred from our modeling results, suggesting that spontaneous subduction of the thicker plate is highly unlikely. Only IBM falls into the OPS domain based on age pairs at onset. There, the initial state of stress being unknown, as both plates edges have been consumed in subduction, it begs the question as to whether the old Pacific plate sunk spontaneously. Now, the question is as follows. To what extent are the rheological parameters and the characteristics of subduction initiation satisfied? Beyond the unrealistic values of crustal densities and brittle properties, the expected sinking rates and asthenosphere rise are high. The slab typically reaches a 200 km length within ∼1 Myr, so the remnants of the resulting “forearc crust” should be restricted to a very short time span. The argument of Arculus et al. (2015) that new findings of juvenile 52–48 Myr old oceanic crust in the Amami-Sankaku Basin that are far away from those already found in the IBM forearc (so-called “forearc basalts”) and the confirmation that the spontaneity of the subduction initiation and wide extent of the asthenosphere invasion was refuted by Keenan and Encarnación (2016) since younger juvenile oceanic crust cannot be used as a test for early uplift in presubduction initiation basement rocks. A second argument comes from the boninitic nature of the primary embryonic arc combining MORB and slab-derived hydrous fluid signatures (Ishizuka et al., 2006). Those boninites erupted between 51 and 44 Ma in the Bonin present arc and forearc (Ishizuka et al., 2011; Reagan et al., 2019) and between 48 and 43 Ma in the Mariana present forearc (Johnson et al., 2014). This time span appears incompatible with the swiftness of the processes required in our models. An alternative geodynamic scenario satisfying both the magmatic and the tectonic constraints has been proposed by Lallemand (2016). The kinematic change of Pacific plate motion following the Izanagi slab break off at approximately 60–55 Ma (Seton et al., 2015) created enabling conditions for convergence across a major TF or fracture zone. Compressive deformation progressively localizes until subduction starts at approximately 52 Ma. At approximately the same time, the occurrence of the Oki-Daitō plume has produced the splitting of the remnant arcs brought by the younger plate (Ishizuka et al., 2013). Oceanic basalts, called FABs (forearc basalts), spread along the axes perpendicular to the nascent trench (Deschamps and Lallemand, 2002, 2003). Boninites erupt as soon as hydrous fluids from the subducting plate metasomatize the shallow asthenosphere beneath the newly formed oceanic crust. Later, tectonic erosion along the IBM trench removes the frontal 200 km of the upper plate, exposing in the present forearc basalts and arc volcanics initially formed far from the trench.

As observed along the Hjort trench, subduction may start as soon as compressive tectonic forces are applied across a TF, but one should note that the subducting plate is the youngest there (Table 1; Abecassis et al., 2016). IBM and Yap cases likely fall in this “forced” category as mentioned above, but we lack direct field evidence, as they were both deeply eroded along their overriding edge.

5.3 Failure of old plate sinking is not excluded

Among the numerous parameter values tested in this study, especially those within reasonable (green) ranges, we have observed that most of them led to incipient subduction of either the young or the old plate but failed soon after (Fig. 6a to e). We have compiled in Table 1 several cases of potential subduction initiation along TFs or fracture zones that either failed (Romanche, Gagua, Barracuda, Tiburon, Saint Paul and Owen) or were just initiated, so we still ignore how it will evolve (Matthew and Hunter, Mussau, Macquarie, and Gorringe). The advantage of studying the aborted cases is that we still have access to the deformation that accompanied subduction initiation, and compression was always recorded in the early stages (see the references in Table 1). These incipient subduction areas are either restraining bends along transform faults or underwent changes in plate kinematics from strike-slip to transpression. A major limiting factor is the cooling of adjacent plates, as the distance from the spreading center or the plume increases, inhibiting their flexural capacities.

6 Conclusions

We perform a large set of 2-D thermomechanical simulations to study the conditions of spontaneous sinking of the older plate at a TF by investigating broad intervals of plate ages and by paying special attention to the mechanical parameter ranges allowing for OPS. OPS is simulated notably if the oceanic crust is dense and mechanically soft far away from the TF on both sides of the plate boundary. Our results confirm that the OP resistance to bending and the YP thickness are the most significant factors preventing OPS. Reducing the brittle properties of the oceanic lithosphere is thus the most efficient way to trigger OPS, compared to a softening by lowering the ductile strength, imposing a hot thermal anomaly or reducing the asthenospheric viscosity. When these extreme conditions are imposed, two processes of OPS are obtained, depending mainly on the assumed YP thickness. They can be summed up as (1) an OP rapid sinking that is decoupled from the YP kinematics and associated with a significant rise in the asthenosphere toward the subducting slab hinge and (2) a dragging of the YP by the sinking OP that is considered a two-sided subduction mode. In all cases, whatever the mode, OPS occurs in less than 1.5 Myr, that is, in an extremely short time span, and only if the initial mechanical setup is adjusted beyond reasonable limits for at least one key thermomechanical parameter. In addition, we find that neither the thermal structure and blurring of the transform fault area nor a plume head impact are able to affect OPS triggering in our modeling setup. Our study highlights the predominant role of a lithospheric weakening to enlarge the combination of plate ages allowing for OPS.

From the parametric study, we conclude that OPS cannot be simulated for a realistic combination of mechanical parameters. By comparing our modeling results to the most likely natural cases where spontaneous subduction at a TF has been invoked, we find that even though extreme mechanical conditions were assumed, the plate age setting at Yap should prevent OPS. Regarding the case of Izu–Bonin–Mariana, in addition to the weakness of geological arguments, the kinetics of subduction initiation, evidenced by geological records, is not compatible with the catastrophic mode systematically simulated in our experiments. We finally conclude that the spontaneous instability of the thick OP at a TF is an unlikely process of subduction initiation in modern Earth conditions.