On earth, the magnitudes of depth-integrated compressive differential stress is in the order of 3–15 TN/m (Cloetingh and Wortel, 1986; Coblentz and Richardson, 1996; Ghosh et al., 2013; Naliboff et al., 2009; Richardson et al., 1979). We emphasise here that this value can vary from one study to another because some used the integral of the maximum horizontal deviatoric stress which is half of the differential stress used in this study ( σ 1 − σ 3 = 2 σ I I ' ;   see Schmalholz et al., 2019; Candioti et al., 2020). At continental passive margins, sources of compressive stresses include: ridge push (e.g., England and Wortel, 1980), gravitational potential energy (GPE, e.g., Pascal and Cloetingh, 2009), tectonic forcing (Bird, 2017), or mantle convection (Ghosh et al., 2013; Kendall and Lithgow-Bertelloni, 2016) (Figure 1). Ridge push arises from lithostatic pressure related to the elevation of the hot mid-ocean ridge above the cooler ocean basins surrounding it (e.g., Forsyth, 1973; Turcotte and Schubert, 2014). Although the contribution of each mechanism is unclear (Swedan, 2015), ridge push represents integrated differential stress values between 1 and 5 TN/m (Mueller and Phillips, 1991; Swedan, 2015; Mahatsente, 2017), with an average value in the order of 3.5 TN/m for a 75 Myr oceanic lithosphere (Mahatsente, 2017). However, in some places, stress arising from the oceanic plate are of equal magnitude as the GPE from thick continental lithosphere (e.g., Tibetan plateau, 7–12 TN/m; Molnar and Lyon-Caen, 1988; Molnar et al., 2015; Schmalholz et al., 2019), suggesting that in addition to ridge push, other forces must contribute to oceanic plate motion (Flesch et al., 2001; Ghosh et al., 2006; Naliboff et al., 2009). The horizontal shear tractions induced by mantle flow (Ghosh et al., 2013), could generate an additional integrated differential stress in the order of 3–6 TN/m (Kendall and Lithgow-Bertelloni, 2016). Consequently, a combination of ridge push and shear traction could generate stress levels in the order of 4–11 TN/m, with an average of ∼8 TN/m for a 75 Myr oceanic lithosphere. Larger integrated forces could only be reached if additional tectonic processes are considered. These are related to horizontal forcing at subduction edges (van Summeren et al., 2012; Bessat et al., 2020), mostly driven by slab pull, which estimates are of the order of 10–20 TN/m (e.g., Fowler et al., 1990; Turcotte and Schubert, 2014). However, as noted in Turcotte and Schubert (2014), the trench pull force is largely balanced by the frictional resistance at the contact between the subducting and overriding plate (Figure 1F).

As loads of 16 TN/m are needed to initiate subduction at passive margins (Zhong and Li, 2020) additional forces are required next to ridge push for the formation of a new subduction zone. It is, therefore, expected that stress levels in excess of ridge push would be relaxed through deformation of the mid-oceanic ridge, where the lithospheric strength is low. The strength of spreading ridges is subject to considerable uncertainty because the magma supply mechanisms are not sufficiently well understood and quantified (Luttrell and Sandwell, 2012). However, inferences from geochemical (Beaussier et al., 2019), geophysical (Luttrell and Sandwell, 2012) and numerical modelling (Husson, 2012) studies suggest that there is a strong correlation between spreading rate and thermal thickness of the lithosphere at mid-ocean ridges. When spreading rates are lower than ca. 1.5 cm yr⁻¹, the melt concentration is particularly low, and the thermal thickness of oceanic lithosphere (isotherm 1,300°C) at the mid-ocean ridges can be as large as 55 km (Husson, 2012). Therefore, an ultra-slow spreading ridge is of similar strength than a 5–8 Myr old oceanic lithosphere. This is in agreement with stress predictions associated with oceanic lithospheric folding suggesting that an ultra-slow mid-oceanic ridge could support horizontal loads of 8 TN/m (e.g., Indian ocean, Gerbault, 2000).

In order to test our hypothesis, we perform numerical modelling simulating shortening of an oceanic basin and the adjacent continental margins. These models do not include melt production processes at mid-ocean ridges but are designed to capture the essential deformation features that characterize the shortening of an oceanic basin of a given age, as a consequence of the interplay between plate cooling and stresses arising from far-field forcing. Ultimately, they serve as a basic model to determine the stress and strain-rate boundary conditions required for subduction initiation at passive margins or within the oceanic domain.

The modelling approach of this study is two-fold (Figures 2A,B). We first run a simplified model setup that adopts a laterally uniform distribution of physical properties within the oceanic lithosphere, where the age of the oceanic lithosphere is everywhere the same (Figure 2A).

In step two, we incorporate the age dependent thickening and lengthening of the oceanic lithosphere from an extinct spreading center to the passive margin (Figure 2B). Hereafter, the extinct spreading center will be referred to as the “ridge”. The results of the parametric study (step 1) facilitate obtaining favourable boundary conditions for subduction initiation at passive margins in step 2, where the role of lithosphere thickness variation which varies through time with cooling was investigated.

The model domain is a section of 3,000 × 500 km and the numerical resolution is 1 × 1 km in both dimensions. The model top boundary is represented by a true free surface (Duretz et al., 2016b). Erosion and sedimentation have been implemented following a kinematic approach (e.g., Candioti et al., 2020) where erosion and sedimentation are implemented above or below a base level fixed at 0 km. In all models, the accommodation space is filled with a sedimentary material composed of calcite, which is a rheologically weak lithology.

Models of step 1 comprise a 2000 km wide oceanic plate flanked by two continental plates on either side (Figure 2A and Table 1). The continental lithosphere consists of an 18 km thick granitic upper crust, a 12 km thick feldspathic middle crust, a 5 km thick granulitic lower crust and a lithospheric mantle with a fixed thermal thickness ( hl ). The oceanic lithosphere entails an 8 km thick crust and a lithospheric mantle with a thermal thickness that depends on its age following the plate cooling model used in Auzemery et al. (2020). The passive margin is characterized by a crust that thins progressively towards the ocean over a distance of 150 km. The thermal base of the lithosphere at the passive margin is defined by a linear interpolation between the ocean and the continent. Deformation is governed by frictional, dislocation, diffusion and Peierls creep equations, with parameters displayed in Table 1. Following previous studies (Zhong and Li, 2019; Candioti et al., 2020; Kiss et al., 2020) we account for thermal softening instead of pre-defined strain softening mechanisms to allow for shear localization leading to subduction initiation. The initial geotherm is computed assuming steady-state conditions and accounts for different radiogenic heat productions in each layer and a constant asthenosphere temperature of 1,330°C.

For each model, we varied the thermal thickness of the oceanic lithosphere and the convergence rate at the boundary, which both influence the stress level within the oceanic lithosphere. The first model setup (uniform thickness) comprises a 2,000 km wide oceanic plate flanked by two continental plates on either side (Figure 2A). 14 ages of oceanic lithosphere were tested, at 10 Myr intervals, counting from 0 to 140 Myr. For each age of oceanic lithosphere, a wide range of convergence rates were tested to arrive at a threshold value for which we observe a change from intra-oceanic subduction to subduction at the passive margin. We carried-out two sets of numerical experiments with two different thermal thicknesses (160 and 180 km) of continental lithosphere fixed in time, representing continental lithosphere of Palaeozoic and Proterozoic age, respectively (Artemieva, 2009).

Based on the results of step one, the thermal thickness of the continental lithosphere has been fixed at 160 km in stage two, as this thickness facilitates subduction initiation at passive margins. The thermal thickness and the length of the oceanic domain depend on the age of the lithosphere along the x-axis, as functions of distance from the former ridge. Therefore, by assuming that the ocean was formed at an average half-spreading rate of 1 cm yr⁻¹ (in accordance with a slow spreading rate, e.g. Dick et al., 2003), the oceanic lithosphere gets 10 Myr older every 100 km (counting from the ridge) and the width of a 30 Myr oceanic basin is 600 km. We carried out three sets of numerical experiments with three different ages of oceanic basin namely 30 and 60 and 90 Myr. For each, the tested parameters include the duration of the cooling period at the start of the experiment and the convergence rate at the boundary during subsequent shortening. The implementation of a cooling period prior to shortening is necessary in order to obtain scenarios where subduction nucleates at the passive margin.

The model top boundary is a true free surface (Duretz et al., 2016b). A constant inward normal velocity ( + V i n ) is applied to the left boundary of the model (Figure 2A) to simulate horizontal tectonic force loading arising from the continent (Duretz et al., 2016a). The right boundary remains fixed. In order to satisfy mass conservation, an outflow velocity ( − V o u t ) is distributed at the base of the model and on the sides from 200 km deep to the base of the box (Figure 2A). The outflow is proportionally distributed over the length of each boundary such as: Σ V i n H i n = Σ V o u t   H o u t

These models with different thermal ages of the oceanic lithosphere have been subject to horizontal velocity boundary conditions. For each age of oceanic lithosphere, a wide range of convergence rates were tested to arrive at a threshold value for which we observe a change from intra-oceanic subduction to subduction at the passive margin.

With the aim of monitoring stress magnitude through time, at every location along the models, we compute for each time step the total integrated stress in the lithosphere I over the thermal thickness of the lithosphere   h l : I = ∫ 0 h l 2 τ I I   d z     = ∫ 0 h l ( σ 1 − σ 3 )   d z     where τ I I   is the square root of the second invariant of the deviatoric stress tensor and ( σ 1 − σ 3 ) is the differential stress (Schmalholz et al., 2019). Though h l corresponds to the 1,330°C isotherm, we use the mechanical base of the lithosphere h m , which is fixed at the 1,000°C geotherm (∼base of the strength envelope, e.g., Figure 2A), for a better visualization of the modelling results.

A stress analysis is performed by computing the integrated differential stress as a function of time (Figure 3) at three positions: the model boundary, the passive margin, and the oceanic lithosphere (Figure 2A). The integrated stress for the oceanic lithosphere is an averaged value over the oceanic domain. At each position, subduction initiation is assumed to start at the time where the integrated stress has reached the highest value and is followed by a sudden stress drop at ts, the time of subduction initiation (Figure 3). Figure 3 is an example of how the integrated stress I varies along the model (for a given set of model parameters). The variation of I describes the state of stress in the lithosphere in terms of increase or reduction of stress relative to the integrated boundary stress at the left side of the model. As such we interpret that when I at the passive margin or in the oceanic domain is lower than at the boundary, the stress is dissipated somewhere, for instance by the presence of ductile shear zones. If I is higher, we interpret that lithosphere undergoes a stress loading until stresses release due to visco-plastic thickening of the plate.

Models were run for a range of convergence velocities, ages of oceanic lithosphere, and thicknesses of continental lithosphere. In each model, the integrated stress at the model boundary was calculated and plotted in various domain diagrams (Figure 4). We also use the evolution of integrated stress through time to explain how it controls the lithosphere dynamics (Figure 5C). A total of 150 models were run in order to establish a robust stress transfer model for the oceanic basin.

In this section, we present the results of the first set of experiments that shows the relationship of convergence rate and the locus of subduction initiation as a function of the age of the oceanic basin (Figures 4A,C). Throughout this study we define t0 as the age of the oceanic lithosphere at the start of the model and ts as the age of the oceanic basin when subduction initiates. For each experiment, we calculate the value of integrated stress I b at the onset of subduction initiation ts (Figures 4B,D). If subduction initiates at the passive margin the value is plotted as a red dot whereas it is shown in blue, when subduction initiates within the oceanic domain.

Figures 4A,C indicate that for young oceanic lithospheres (age < 50 Myr) subduction initiation at a passive margin is only feasible at low convergence rates. The upper bound value of convergence rate for subduction at passive margin (max v ) is relatively constant (∼ v = 1.5 cm yr⁻¹, Figure 4A) and thus independent of the age of the oceanic lithosphere. Moreover, this limit is relatively lower for 180 km thick continental lithosphere (∼ v = 0.9 cm yr⁻¹, Figure 4C).

In contrast, for oceanic lithosphere older than 50–60 Myr, the upper bound values of convergence rate for subduction initiation at passive margin varies with the age of oceanic lithosphere at the start of the model as well as the thickness of the continental lithosphere. Intra-oceanic subduction is consistently observed for higher convergence rates across all tested ages of oceanic lithosphere (Figures 4A,C). Moreover, the range of velocities leading to subduction initiation at a passive margin is much larger in case of a 160 km thick continental lithosphere (Figure 4A), than for a 180 km thick continental lithosphere (Figure 4C), suggesting that subduction at passive margin is less likely in case of thick continental lithospheres. It shows again that for similar kinematic conditions, strain localization at passive margins is controlled by the strength of the continental lithosphere, regulated through its thermal thickness.

From the distribution of modelling results, we delineate upper and lower bounds of stress levels for subduction initiation (max. and min stress, Figures 4B,C). The results show that, overall, subduction at passive margins requires less stress than intra-oceanic subduction suggesting that subduction at passive margins is possible for stress levels lower than that within the oceanic plate (Figure 4B, blue dash-dot line), which is predicted for oceanic lithospheres older than ∼40 Myr (Figure 4B). Consequently, with cooling of the oceanic lithosphere, intra-oceanic subduction becomes particularly intricate and requires large amounts of integrated stress.

The minimum integrated differential stress necessary for subduction nucleation at a passive margin (red dashed line at 60 Myr, Figure 4B) is in the order of 30 TN/m, which corresponds to an age of 40–50 Myr for the oceanic lithosphere at the time of subduction. The stress magnitude required for subduction at passive margins is sensitive to the thermal thickness hl of the continental lithosphere. The minimum stress for subduction initiation at passive margins ranges from 30 to 36 TN/m for hl = 160 km and from 40 to 50 TN/m for hl = 180 km (Figures 4B,D), which corresponds to age limits of oceanic lithosphere of 40 Myr for hl = 160 km and 60 Myr for hl = 180 km. This result shows that among all the parameters analysed in this study, the thermal thickness of the continental lithosphere is probably most important. For similar ages of subduction initiation at passive margins, the disparity in stress is due to a variation in convergence rate, with low convergence rate models requiring lower stress levels for subduction initiation at passive margins (red dashed line, Figure 4B).

Subduction initiation at passive margins requires high loading level of shear stress acting on the margin but also low levels of shear stress in the oceanic lithosphere. Our modelling results show that convergence velocity is an important parameter that regulates stress in the lithospheric layers and therefore controls the locus of deformation. Models with low convergence rate (∼0.8 cm yr⁻¹) predict strain localization and subsequent subduction initiation at the passive margin (Figure 5A, 6.6%BS), whereas models with high convergence rate predict intra-oceanic subduction (Figure 5B, 7%BS). These results show that stress loading preferentially occurs at passive margins for cases of low convergence rate (Figure 5C), because deformation is distributed over many small structures within the brittle layer that accommodate low amounts of strain and is even more distributed within the ductile layers of the models. Consequently, distributed deformation (black curve) leads to low shear rates and low heat production (Figure 5A). The heat produced by viscoplastic dissipation ( Hs , Figure 5A) is then efficiently diffused within the viscous layer (red curve, Figure 5A, 6.6%BS). Therefore, homogeneous distribution of deformation and shear heating at low strain-rate limits the magnitude of stress in the oceanic domain and thus prevents stress loading and failure during cooling of the lithosphere (Figure 5C). When the oceanic lithosphere reaches a certain thickness (after t = 30–35 Myr), it barely deforms, but acts as a buttress and shortening leads to deformation of the passive margin until subduction initiates (Figure 5A, 8.5%BS).

The passive margin lithosphere consists of brittle layers in the crust and the underlying lithospheric mantle where deformation localises in shear bands, which eventually link-up to form a through-going shear structure (Figure 5A). This moment of formation of the incipient subduction plate boundary allowing for underthrusting of the oceanic plate correlates with a significant stress drop as shown in Figure 5C.

In contrast, for models with relatively high convergence rate deformation is localized in the oceanic domain (Figure 5B, between -500-0 km). The increase of integrated stress is largest within oceanic lithosphere for both young and old oceanic basins (Figure 5D, t = 20 Myr). Higher strain rates within shear band type-structures lead to localised shear heating ( Hs ) triggering shear-localization by thermal softening in the ductile layer and the formation of a subduction plate boundary (Figure 5B).

In this section we present the results for three sets of experiments defined by the age of the oceanic basin (30 Myr, 60 Myr and 90 Myr; Figures 6A–C) that include age dependent thickness variations within the oceanic lithosphere (Figure 2B). Figure 6 delineates the mode of subduction as function of the cooling period implemented at the start of the experiments and the convergence rate during the subsequent shortening period. Overall, these modelling results show that a cooling period of at least 5 Myr before the start of shortening (Figure 6A) is necessary to prevent underthrusting at the position where the thermal thickness of the oceanic lithosphere is smallest, representing in simplified forma former ridge. Consistent with results derived from step 1, we infer that subduction initiation at a passive margin is only feasible for low convergence rates ( v < 0.9 cm yr⁻¹, Figure 6A). However, different to the models with a uniform oceanic lithosphere, Figure 6 shows that subduction initiation is more favourable in case of a young oceanic basin rather than an old one (Figures 6A vs Figures 6B,C). To understand this result, we present in Figure 7 two modelling results with similar boundary conditions (tc = 5 Myr, v = 0.3 cm yr⁻¹) but with different ages of the oceanic basin (age = 30 and 60 Myr; Figures 6A,B) at the onset of shortening.

For both models, the oceanic lithosphere shows an increase in thickness and integrated stress from the basin centre to the continental lithosphere (Figures 7A,B), which is largest for the older (60 Myr) oceanic basin (Figures 7A,B). In case of a 30 Myr oceanic basin, the level of integrated stress at the ridge is close to the value at the margin (∼16–18 TN/m, black curve in Figure 7A). Shortening of such young lithosphere at a slow rate of 0.3 cm yr⁻¹ leads to strain localization at the former ridge early in the deformation history (Figures 7A, 7 Myr). However, in the later stage (Figure 7A, 30 Myr), strain is rather distributed, because lateral thermal thickness variations get reduced as the oceanic lithosphere cools faster at the ridge compared to the passive margin. Under such conditions, the subsequent development of a shear-zone at the base of the continental crust leads to underthrusting of the oceanic plate under the margin (Fig. 7a, 46 Myr). In comparison, in case of a 60 Myr old oceanic basin the level of integrated stress at the ridge is significantly lower than at the passive margin ( ΔI = 10 TN, Figures 7B, 7 Myr).

Consequently, strain localises at the mid-oceanic boundary which is much thinner (Fig. 7b, 45 Myr), resulting in intra-oceanic under-thrusting (Fig. 7b, 45 Myr). These results suggest that the strength contrast (expressed in this section by the thermal thickness) between the lithospheres at the centre of the oceanic basin and the passive margin controls the locus of subduction initiation. Therefore, a young oceanic basin with minor variation in strength between the extinct ridge and the margin requires only 5 Myr of cooling prior to convergence to permit subduction at a passive margin. In contrast, an old oceanic basin would require a significant period of cooling (tc > 20 Myr) to reduce the strength differences between the lithosphere at the ridge and at the passive margin to allow for subduction initiation at the passive margin (Figure 6B).

We have set to investigate what conditions and mechanisms lead to intra-oceanic subduction vs passive margin subduction (Figure 8). Because mid-oceanic ridges have a thinner and weaker lithosphere compared to passive margins, they represent a priori preferential locations for subduction initiation (e.g. Maffione et al., 2015; Agard et al., 2016). Moreover, the large difference between ridge push force (∼1–5 TN) and stresses arising from tectonics, or topography related gravitational potential energy (5–10 TN) suggests that subduction should predominantly initiate at mid-oceanic ridges, which is consistent with our modelling predictions for fast convergence rates. Additionally, this modelling study suggests that subduction initiation at passive margin is physically feasible under restricted conditions.

First, subduction initiation at a passive margin requires processes that regulate stress levels in the oceanic lithosphere. Our modelling results show that rheological conditions and processes favouring distributed deformation within the ductile part of the oceanic lithosphere is key for reducing stress levels and de-localising deformation (Figure 8B). As such localised deformation within the brittle layer, expressed as shear band type structures, fail to propagate into and through the ductile layer to form a subduction plate boundary. This behaviour is tied to low convergence rates, i.e., below 0.9 cm yr-1 (Figure 8B). Similar results have been obtained by Gülcher et al. (2019) and Qing et al. (2021), who observed in their numerical models, reactivation of multiple detachment faults upon their inversion but did not produce intra-oceanic subduction. When deformation does not localized in oceanic domain, the vertical rheological decoupling at the margin allows for the development of a long-lasting shear zone where stresses are relaxed through the formation of a decollement that propagate through the mantle. Our results predict that even young (5 Myr old) oceanic lithospheres can support tectonic stress of up to 12 TN/m underlining its long-term stability. As such, subduction initiation along spreading ridges is largely favored by warm ridge and fast convergence rate (Figure 8C, see also Qing et al., 2021) and/or the implementation of lithosphere-scale pre-existing weak zones (Maffione et al., 2015) to localize deformation in the mantle lithosphere.

Second, subduction initiation at a passive margin also entails the transfer of deformation from the ridge to the margin during the cooling of the oceanic lithosphere (Figure 7A). Previous studies have emphasized that deformation at the passive margin through spontaneous margin collapse is an unlikely mechanism for subduction initiation, because stresses acting on the passive margin lithosphere are never at yield (Cloetingh et al., 1984; Mueller and Phillips, 1991). It is thus more likely that the nucleation of a subduction zone at a passive margin occurs upon additional external forcing to reach stress levels of at least 16 TN/m (Zhong and Li, 2020) to induce failure of the passive margin lithosphere. Numerical modelling studies simulating the shortening of an oceanic basin (Auzemery et al., 2020; Candioti et al., 2020; McCarthy et al., 2020; Zhong and Li, 2020) emphasise that such stress levels can only be supported by thick oceanic lithospheres, suggesting that subduction initiation is only feasible in case of shortening of an old oceanic basin. Therefore, subduction at a passive margin could only happen after a relatively long period of cooling at an extinct ridge prior to plate convergence (>60 Myr, Candioti et al., 2020; McCarthy et al., 2020). In comparison, our models simulating the effects of an extinct and cooling mid-oceanic ridge suggest that subduction initiation at passive margins is also possible in case of slow shortening of a young oceanic basin (<30 Myr), following a short period of cooling (>5 Myr). Although our approach assumes an initial short period (5 Myr) of cooling where no external forcing (convergence) is applied, the results predict a more reasonable time scale for the development of a subduction zone at a passive margin, which amounts to 45 Myr from the moment that shortening is applied.

Third, subduction initiation at continental margins does not only depend on the age of the oceanic lithosphere (e.g. Auzemery et al., 2020; Zhong and Li, 2020) but also on the thermo-tectonic age of the rifted continental lithosphere (Figure 4, see also Nikolaeva et al., 2010). We show that subduction initiation is favoured for thermally thin (Palaeozoic or younger) continental lithospheres (<160 km) over cratons (>180 km) (Artemieva, 2009). This could explain why subduction does not exist along the north and south Atlantic margin, where the continental lithosphere is particularly thick (Tesauro et al., 2013) and might also be the reason for the proposed transference of subduction from the paleo-Tethys, to the neo-Tethys, along the margin of the thinned Cimmerian micro-continent (Wan et al., 2019).

Similar to Candioti et al. (2020) and Kiss et al., 2020, the stress levels obtained in our study are larger than plate boundary forces and represent upper bonds. Thus, other pre-existing lithospheric structures or weakening mechanism are generally suggested to explain subduction initiation at a passive margin (e.g., Stern and Gerya, 2018).

First we note that, the level of integrated stress required for subduction initiation depends mainly on the strength of the continental lithosphere that depends largely on its composition, thermal regime and the presence of weakening mechanisms (Cloetingh et al., 2005) and fluids (Regenauer-Lieb et al., 2001). In our study, the stress limit also varies with rheological layering and thermal state at the margin (Auzemery et al., 2021) and the earlier mentioned 20 TN/m limit is only valid for a 4 layer-160 km thick-continental lithosphere. This explains differences with similar studies (e.g., Zhong and Li, 2019), in which the crust is weaker and thermal thickness at the margin is significantly smaller than in our models.

Second, in our models, the sensitivity of deformation pattern to the choice of softening parameters is too high. Therefore, to not prescribed the locus of deformation, the magmatic (Gerya and Meilick, 2011), hydro-mechanical (Dymkova and Gerya, 2013; Schmalholz et al., 2020) and grain boundary weakening (Bercovici and Ricard, 2014) were not implemented. Parameterized strain softening was also not used because it is not a transient mechanism and the strain-dependent weakening limits are not well calibrated (Jaquet and Schmalholz, 2018). In addition, our models do not include inherited structures such as detachment faults that are prone to reactivation at stress levels (Maffione et al., 2015) that are well below those required for the development of new faults. Furthermore, a very thin oceanic lithosphere is used to simulate an extinct mid-ocean ridge. As such, melt- and fluid-related processes at mid-oceanic ridges, which are yet poorly constrained because of scars high resolution observations and rheology data from natural systems (Bickert et al., 2021) are not taken into account. These have been discussed in models deploying lithosphere extension (e.g. Ligi et al., 2008; Dannberg et al., 2019), which is different to our setup, which focuses on contraction of an oceanic ridge and its continental margins. The mechanisms discussed above could greatly reduce the stress limits after deformation localize at a passive margin or at a mid-oceanic ridge. However, the doubts and uncertainties associated with the choice of the weakening parameters would affect too much the results of the parametric study and a less complex approach was adopted.

Lastly, the discrepancy between analytical and numerical models could be explained by the very high uncertainties about the levels of integrated stress present within the lithosphere. Traditionally, the magnitudes of integrated stress used to explain plate tectonic processes are compared with analytical calculations made on plate boundaries (subduction, ridge-push; e.g., England and Wortel, 1980). However, these models are extremely unclear as they ignore multi-scale and multi-physics processes such as magmatic accretion or mantle dynamics (e.g., Husson et al., 2015). Besides, in nature, additional stress sources are necessary in order to explain regional-to local-scale stress (e.g. Heidbach et al., 2007) and deformation patterns (Gerbault, 2000; Cloetingh et al., 2005). Furthermore, it is common that forces are calculated using the integral of the horizontal deviatoric stresses (Ghosh et al., 2006; Mahatsente and Coblentz, 2015) that are lower by a factor of two (see Schmalholz et al., 2019) than the differential stress usually used to calculate lithosphere forces.

Our modelling results together with previous modelling studies (Cloetingh et al., 1989; Gurnis et al., 2004; Hall, 2019; Nikolaeva et al., 2011; Stern and Gerya, 2018) suggest a mechanism that favors forced subduction initiation at passive margins. In particular, we find that slow convergence (<0.9 cm yr⁻¹) over a long period of time (∼40 Myr) is critical in this context. Forces at play may include ridge push (e.g., Forsyth, 1973; Vlaar and Wortel, 1976) GPE from adjacent high areas (Ghosh et al., 2006; Marques et al., 2013; Pascal and Cloetingh, 2009), transference from an existing subduction zone (Baes et al., 2018; Duarte et al., 2013) or mantle flow (Candioti et al., 2020 and references therein). Associated with this far-field tectonics (Figure 8A), the partial starvation of basaltic melt underneath an ultraslow spreading ridge can lead to the gradual increase in thermal thickness of the oceanic lithosphere by 40–50 km (Husson, 2012) resulting in a period where stresses transmitted from the extinct ridge affect the passive margin (Figure 7).

In the frame of the above, the Alpine Tethys is a particular good analogue where the geological record demonstrates ultra-slow spreading (Lagabrielle and Cannat, 1990) and subduction initiated at the passive margin (Manzotti et al., 2014b; Marroni et al., 2017; Schmid et al., 2004). The latter process was probably associated with low convergence rates affecting the Adriatic crust of Palaeozoic (∼320 Ma) tectono-thermal age (Castellarin and Cantelli, 2010). The actual age of subduction initiation in the Alps is debated, but high pressure rocks found in the Adriatic margin constrain an age for underthrusting of ca 90–80 Ma (Manzotti et al., 2014a and references therein). However, several authors suggest an earlier subduction initiation at ca. 130 Ma, related to the change in Africa-Europe convergence controlled by the opening of the south Atlantic Ocean (Handy et al., 2010; van Hinsbergen et al., 2020; Le Breton et al., 2021). Based on our modelling results, we argue that subduction initiation in the Alps was a long-lasting (40–50 Myr) process that required a long period of very slow convergence allowing for strengthening of the ridge and the stress build-up at the ocean-continent transition. We note that the latter is achieved for times when Africa did not move head-on relative to Europe. During this period from 130 to 85 Ma, referred as “Cretaceous Quiet Zone” by van Hinsbergen et al. (2020), plate reorganisation leads to a quasi-absence of convergence between Africa and Europe (c.a. 0.6 cm yr⁻¹ based on Figure 9). We suggest that such configuration would favour a reduction of the spreading rate, the strengthening of the mid-oceanic ridge as well as strain accumulation at the margin, until subduction finally initiates due to an increase in plate convergence rate 90 Ma (Capitanio and Goes, 2006). This scenario implies that at the onset of slow convergence at 130 Ma, the Piemonte-Liguria ocean is 30–40 Myr old, which falls within the range of “favourable” conditions for SI at passive margins (Figure 6A). We furthermore remark that exhumed and serpentinised mantle lithosphere of the Alpine Tethys provided additional favourable mechanical conditions for subduction initiation at stress levels that are compatible with plate tectonic forces (e.g. Candioti et al., 2020).