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 ﬁrst seek candidates from recent subduction initiation events at an oceanic TF that could fulﬁll 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 rhe-ologies, 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 astheno-sphere may rise up to the surface above the sinking old plate, provided that the younger plate remains motionless (veriﬁed for ages ≥ 5 Myr, mode 1). For lower younger plate ages (typically ≤ 2 Myr


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 Published by Copernicus Publications on behalf of the European Geosciences Union.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.
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 un-derwent 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 andLallemand, 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.

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(Zhu et al., , 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 −1 ).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.

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 Solid Earth, 11, 37-62, 2020 www.solid-earth.net/11/37/2020/that may be involved in the process of subduction initiation (Fig. 2).Density (ρ) is assumed to be temperature and composition dependent: 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, ρ ref m is fixed to 3300 kg m −3 , while ρ ref for the oceanic crust and the weak material is varied from one experiment to another (Sect.2.4).

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: where C 0 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): ε = εref (τ/τ y ) n p , where ε is the second invariant of the deviatoric strain rate tensor, ε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 τ = ν ε, the plastic viscosity, ν b , is written as follows: A dislocation creep rheology is simulated using a non-Newtonian viscosity ν d , defined by where B 0 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 2.6 × 10 −21 s −1 , but no maximum cutoff is imposed.
www.solid-earth.net/11/37/2020/Solid Earth, 11, 37-62, 2020  (Ribe and Christensen, 1994;Arcay, 2012).The red dotted line represents the hot thermal anomaly ( T = +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.

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 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 erf −1 (0.9) ∼ 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.64www.solid-earth.net/11/37/2020/Solid Earth, 11, 37-62, 2020 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.

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.

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 c a (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 c a , γ 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 γ TF = 5 × 10 −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, γ TF = γ c = 0.05 and L w = 0 km: the TF and both plate surfaces are made of gabbroic oceanic crust (Table 3).Note that when γ c = γ TF = 5 × 10 −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 c a 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).
Solid Earth, 11, 37-62, 2020 www.solid-earth.net/11/37/2020/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 m a .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 s = τ/σ n , where τ is the shear stress along the fault (Turcotte and Schubert, 1982): where λ is the pore fluid pressure coefficient, ρ w is the water density, and p w is the pore fluid pressure, assuming that p w = ρ w gz 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).
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).

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 −3 ; Bousquet et al., 1997;Tetreault and Buiter, 2014).Crust density in the blueschist facies reaches 3160 kg m −3 , 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 www.solid-earth.net/11/37/2020/Solid Earth, 11, 37-62, 2020 is here a function of temperature through the coefficient of thermal expansion (Table 2).

Activation energy for the crust
The most realistic interval for the crustal activation energy E c a can be defined from experimental estimates E exp a for an oceanic crust composition.Nonetheless, E exp a 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 c a 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 2/(n + 1) (Dumoulin et al., 1999).As sublithospheric small-scale convection yields strain rates by the same order of plate tectonics (∼ 10 −14 s −1 ; 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: where n e is the experimentally defined power law exponent.The activation energy E exp a in the dislocation creep regime is encompassed between the one for a microgabbro, 497 kJ mol −1 (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 −1 (Mackwell et al., 1998, with n e = 4.7 ± 0.6).For a basalt, E (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 exp a = 154 kJ mol −1 , 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 −3 , Fig. 3e).

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.

Numerical code and resolution
The models are performed using the thermo-chemomechanical code of convection developed by Christensen (1992), which is based on an Eulerian and spline finite element method.Conservation equations are solved to obtain Solid Earth, 11, 37-62, 2020 www.solid-earth.net/11/37/2020/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 2 ), 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).

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.

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 −1 ) 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).

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 −1 in simulation S1a, 0 vs. 2, up to ∼ 180 cm yr −1 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).

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).

Transform fault and oceanic crust densities
Densities strongly affect the evolution of the TF system.If the TF weak medium is buoyant (ρ TF = ρ c = 2920 kg m −3 ), 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 (2 < A y < 20 Myr).On the other hand, a heavy material filling the TF gouge (ρ TF = 3300 kg m −3 ) 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).

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.

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 (γ c = 5 × 10 −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 −3 and 5 × 10 −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.

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 γ crit m below which OPS occurs.This threshold depends on the OP age: γ crit m ∼ 0.06 for the plate age pair 10 vs. 40 (simulations S37c to f) but the threshold is higher (γ crit m ≥ 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 γ crit m 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 (γ crit m 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 m a (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 m a 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).
Solid Earth, 11, 37-62, 2020 www.solid-earth.net/11/37/2020/In panel (f), the boundary between the "no subduction" and "OPS" domains corresponds to the relationship A o /A 2.5 y 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.

Ductile strength decrease
The main effect of imposing a decrease in the crust and TF ductile strength (lowering E c a to 185 kJ mol −1 ) 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 c a 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 y , A o ) = 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 c a allows for OPS in this unique case.The YP destabilization and dripping result from the high crust density (ρ c =3300 kg m −3 ) assumed in experiments performed with a reduced E c a .Indeed, in experiments using a usual crustal density (ρ c = 2920 kg m −3 ), YP vertical subduction is obtained instead (simulation S15g, 2 vs. 10, for instance).

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 c a : 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 −3 ).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 eiwww.solid-earth.net/11/37/2020/Solid Earth, 11, 37-62, 2020 ther 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.7Additional 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 5 × 10 −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.

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.

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 −1 , 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 −1 ).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 2 to 10 3 ; 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 2 .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.
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.
If parameter increases (↑) . . . . . .does the force promote (+)   or inhibit 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 o /A 2.5 y 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.4on 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.

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 c a 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.

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.
Solid Earth, 11, 37-62, 2020 www.solid-earth.net/11/37/2020/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.

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. 6fg).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).

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 plumeinduced 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 −1 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.

Weakening of the oceanic mantle lithosphere
Our results show that OPS is strongly facilitated if the lithospheric mantle is weak.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.Grainsize 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 −21 m 2 ) 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 −19 and 10 −16 m 2 , which would hamper high pore fluid pressures and, eventually, plate bending.

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 ex-tremely young (< 7 Myr), the oceanic crust should be very dense (ρ c = 3300 kg m −3 ), as well as drastically weakened (γ c = 5×10 −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 −1 .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 Solid Earth, 11, 37-62, 2020 www.solid-earth.net/11/37/2020/the younger plate (Ishizuka et al., 2013).Oceanic basalts, called FABs (forearc basalts), spread along the axes perpendicular to the nascent trench (Deschamps andLallemand, 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.

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.

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.Marguerite Godard, Martin Patriat, and Andrea Tommasi.We are grateful to Stéphane Arnal, Fabrice Grosbeau, and Josiane Tack for the maintenance and development of the lab cluster of computing nodes on which all numerical experiments were performed.We are also grateful to Anne Delplanque, who drew Fig. 1.We warmly thank Ben Maunder and an anonymous reviewer for their thorough and very constructive comments that significantly improved the article, as well as Susanne Buiter for handling the editing.Review statement.This paper was edited by Susanne Buiter and reviewed by Ben Maunder and one anonymous referee.

Figure 1 .
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.

Figure 3 .
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 ca , 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).
exp a has been recently estimated to 456 kJ mol −1 (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 exp a = 212 kJ mol −1 , n e = 2.4; Ranalli, 1995), and are used to define the lower bound of the "yellow" range for E c a (Fig. 3e).A few experiments have shown that E exp a can be as low as 132 kJ mol −1

Figure 4 .
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 .
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.

Financial support .
This research has been supported by CNRS-INSU (National Institute of Universe Science) program "TelluS-SYSTER" (2015 and 2016).

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

Table 2 .
Constant names and values.

Table 3 .
List of simulations quoted in the text.See the data in the Supplement for the complete simulation list.If one value only is indicated, the oceanic crust (3) in Fig.2is assumed to have the same brittle parameter as the weak material forming domains (1) and (2).bIfonly one value indicated, then ρ c = ρ TF .cWhenone value only is indicated, L w is identical on both plates.dTheweak 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 19 Pa s −1 ).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.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 a Demouchy et al., 2013)e 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).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; (Garel et al., 2014)979)ticity(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) fluidwww.solid-earth.net/11/37/2020/Solid Earth, 11, 37-62, 2020 related weakening.