Here we focus on SZI at magma-poor passive margins that bound embryonic oceans. We define embryonic oceans as marine basins that are small in extent (ca. ≤ 600 km) and mainly floored by exhumed mantle and hyperextended continental crust rather than mature oceanic crust (Chenin et al., 2017). Ancient and recent examples of embryonic oceans are, for example, the Piemonte–Liguria ocean (e.g., McCarthy et al., 2020), the present day Ligurian Basin (e.g., Dannowski et al., 2020) and Red Sea (e.g., Chenin et al., 2017). These ocean basins result from slow-to-ultra-slow spreading rift systems, and the architecture and extent of these oceans is similar (Chenin et al., 2017). Mesozoic rifting separated Europe from Adria and terminated at ca. 145 Ma (e.g., Le Breton et al., 2021). This rifting generated the Piemonte-Liguria ocean basin and the two magma-poor passive margins of Europe and Adria (e.g., Handy et al., 2010; McCarthy et al., 2020; Le Breton et al., 2021). The European and Adriatic passive margins were likely hyperextended and bounded the embryonic Piemonte–Liguria ocean (see Figure 1A). At ca. 85–80 Ma subduction was presumably initiated within the Adriatic continental margin leading to subduction of the Sesia–Dent Blanche unit (e.g., Manzotti et al., 2014, see Figure 1B). Remnants of the Piemonte–Liguria ocean have been incorpated in the Alpine orogenic wedge (e.g., Schmid et al., 2017, see Figure 1C). Whether the European slab is detached or still attached to the European plate is currently disputed (e.g., Kästle et al., 2019).

In this study we build upon our previous work presented in Candioti et al. (2021). Candioti et al. (2021) presented two-dimensional (2D) petrological-thermomechanical numerical models that simulate several, subsequent geodynamic processes in a single and continuous simulation, including thermal convection in the upper mantle. They modelled 1) the formation of hyperextended passive continental margins bounding an embryonic ocean, 2) SZI via thermal softening during convergence of the ocean-margin system, 3) subduction and closure of an embryonic ocean, 4) formation of a collisional orogen and 5) SD of subducted lithosphere. This Wilson-type cycle of embryonic ocean basins comprised 1) a 50 Myr extension period, 2) a 60 Myr cooling period between the extension and convergence phase and 3) a ca. 70 Myr convergence period. Candioti et al. (2021) applied this cycle to the opening and closure of the Piemonte–Liguria ocean and the subsequent formation of the European western Alps (e.g., McCarthy et al., 2020). They mainly focused on the competition between buoyancy and shear forces during the collisional stage of the orogen. When shear forces dominate the collisional stage of the orogen, wedge formation without deep subduction of continental crustal material occured. When buoyancy forces were equally important as shear forces, Candioti et al. (2021) observed deep subduction and exhumation of continental units. The mechanical boundary velocities and the duration of the deformation periods they applied were motivated by plate motion reconstructions from the European and Adriatic plate during the past ca. 180 Myr (e.g., Le Breton et al., 2021). Complementary to their results, we here focus on the subduction stage. Especially, we aim at quantifying the minimum force required for SZI and the role of lithospheric strength in the kind of models presented by Candioti et al. (2021). Note, that we refer to time as the simulated physical time and its unit is Myr. Whenever we refer to actual geological time, we use the unit Ma. Figures 2A-D shows the initial configuration used by Candioti et al. (2021) to model this type of Wilson cycle. A detailed justification of the initial conditions, a full description of the entire model evolution as well as in depth information on the implementation of boundary conditions, petrological-thermomechanical coupling and surface processes (i.e., sedimentation and erosion) is given in the open-access article of Candioti et al. (2021). The models we present here start from the configuration predicted at 110 Myr (end of cooling period) in model history by the reference model of Candioti et al. (2021) (see Figures 2E-H). Similar to the study of Candioti et al. (2021), we parameterize a serpentinisation front propagating across the top 3 km of exhumed mantle material in the embryonic oceanic basin before the onset of convergence. Because of ultra-slow spreading rates, no significant mid-ocean ridge with active magmatism is formed. Therefore, embryonic oceans often lack a mature Penrose-type oceanic crust (e.g., McCarthy et al., 2020). Instead, the subcontinental mantle is exhumed to the surface and in contact with the sea-water. Thus, serpentinisation is likely significant in embryonic ocean basins and presumably plays an important role during subduction of embryonic oceans (e.g., McCarthy et al., 2020). In our models, convergence is induced by applying a constant absolute boundary horizontal velocity of 1.5 cm yr⁻¹. To quantify the minimum driving force and assess the strength of embryonic oceanic lithosphere, we employ a parameterized stress limiter in the mantle lithosphere and the lower crust. In addition to the applied pressure-sensitive Drucker-Prager plastic stress limiter function, this stress limiter is independent on pressure and is used here as a controlling parameter for the plate strength. We systematically tested stress limits of 10, 50, 100, 150, 200, 300, 500 MPa, and 500 GPa. Employing a value of 500 GPa deactivates the stress limiter in the reference model, named NOSLIM hereafter. In addition to the reference model, we here report results for models employing stress limits of 100, 200 and 300 MPa, named SLIM100, SLIM200 and SLIM300 hereafter. These models predict end-member dynamics for subduction initiation and the evolving subduction zone. To test the impact of shear heating and associated thermal softening on SZI we also performed the two simulations SLIM300 and NOSLIM without shear heating.

We consider incompressible materials that slowly creep (no inertia) under gravity and boundary tractions. The applied numerical algorithm solves the continuity equation, the equation for the conservation of momentum and the heat transfer equation in a two-dimensional Cartesian coordinate system, which are given by: ∂ v i ∂ x i = 0 (1) ∂ σ i j ∂ x j = − ρ g i (2) ρ c P D T D t = ∂ ∂ x i k ∂ T ∂ x i + H A + H D + H R , (3) where i and j are spatial indices and repeated indices are summed, v is velocity and x is the spatial coordinate, σ is the total stress tensor, ρ is density and g = [0; −9.81] is a vector including the gravitational acceleration. The density is a function of temperature, T, pressure, P, and a given bulk rock composition, C, which is precomputed using a thermodynamic software package (Perple_X Connolly, 2005, 2009). We, thus, include large density variations resulting from mineral phase transitions. For example, at a lithostatic pressure and a temperature corresponding to a depth of ca. 410 km, the olivine–wadsleyite phase transition results in a density increase of ca. 150 kg m⁻³ (see Candioti et al., 2020, Figures 2B,E). Such an increase might significantly impact on the subduction dynamics. A detailed description of the implementation of petrological-thermomechanically coupled processes and the applied chemical composition of the model units is given in the two open-access publications Candioti et al. (2020) and Candioti et al. (2021). In Eq. 3, c _P is heat capacity at constant P, D / D t is the material time derivative, k is thermal conductivity and H _R comprise source terms resulting from radiogenic heat production. We include source terms resulting from adiabatic processes (Extended Boussinesq Approximation, e.g., Candioti et al., 2020, for detailed justification), H _A = Tαv _y gρ, where α is the coefficient of thermal expansion and v _y is the vertical component of velocity. Further, we consider contributions from the conversion of rate-dependent and rate-independent mechanical processes into heat, H D = τ i j ε ̇ i j − ε ̇ i j ela , where ε ̇ is the total deviatoric strain rate tensor and ε ̇ ela denotes the strain rate resulting from elastic deformation. The total stress tensor is decomposed into a pressure and a deviatoric part as σ i j = − P δ i j + 2 η eff ε ̇ i j eff (4) δ i j = 0 , i ≠ j 1 , i = j , (5) where δ _ij is the Kronecker-Delta, η ^eff is the effective viscosity, ε ̇ eff is the effective deviatoric strain rate tensor, ε ̇ i j eff = ε ̇ i j + τ i j o 2 G Δ t , (6) where G is shear modulus, Δt is the time step, τ ^o are the deviatoric stress tensor components of the preceding time step. We consider visco-elasto-plastic deformation and additively decompose (so-called Maxwell model) the total deviatoric strain rate tensor into contributions from rate-dependent plastic deformation (i.e., effectively viscous behaviour by dislocation, diffusion and Peierls creep), elastic and rate-independent (brittle) plastic deformation as ε ̇ i j = ε ̇ i j ela + ε ̇ i j pla + ε ̇ i j dis + ε ̇ i j dif + ε ̇ i j Pei . (7)

The viscosity for dislocation, diffusion and Peierls creep are a function of the respective second strain rate tensor invariant, ε ̇ II dis, dif,Pei = 1 2 ε ̇ i j dis, dif,Pei 2 . The dislocation and diffusion creep viscosity are computed as η dis, dif = 2 1 − n n 3 1 + n 2 n A − 1 n ε ̇ II dis, dif 1 n − 1 d m exp Q + P V n R T f H 2 O − r n , (8) where d is grain size and the ratio in front of the prefactor A results from the conversion of experimentally derived flow laws into a tensor formulation of these flow laws (e.g., Schmalholz and Fletcher, 2011). The quantities A, n, m, Q, V, f H 2 O and r are experimentally determined material parameters which are characteristic for each creep law (see Table 1). To increase the readability we have omitted the superscripts “dif/dis” for these parameters. Note that for diffusion creep the n = 1 and for dislocation creep m = 0. Hence, dislocation and diffusion creep are essentially insensitive to grain size and strain rate, respectively. We further include the experimentally derived Peierls flow law by Goetze and Evans (1979) expressed in the regularized form of Kameyama et al. (1999) and calculate the Peierls viscosity as η Pei = 2 1 − s s 3 1 + s 2 s A ̂ ε ̇ II Pei 1 s − 1 , (9) employing an effective stress exponent s = 2 γ Q R T 1 − γ . (10) A ̂ in Eq. 9 is calculated as A ̂ = A Pei exp − Q 1 − γ 2 R T − 1 s γ σ Pei , (11) where A ^Pei, γ and σ ^Pei are material parameters. We control rate-independent (brittle) plastic deformation by a deviatoric stress limiter function F = τ II − τ eq (12) τ eq = P sin ⁡ φ + C cos ⁡ φ , P sin ⁡ φ + C cos ⁡ φ < S S , P sin ⁡ φ + C cos ⁡ φ ≥ S , (13) where, τ II = 1 2 τ i j 2 is the second invariant of the deviatoric stress tensor, φ is the internal angle of friction, C is the cohesion and S is a stress limiter which is a parameter that is not experimentally derived but a variable used to control the material strength. The plastic viscosity at the equivalent deviatoric stress, τ ^eq, is computed as η vep = τ eq 2 ε ̇ II eff , (14) where ε ̇ II eff = 1 2 ε ̇ i j eff ε ̇ i j eff . In Eq. 4, the effective viscosity is either the quasi-harmonic average of the visco-elastic contributions η eff = 1 G Δ t + 1 η dis + 1 η dif + 1 η Pei − 1 , F < 0 η vep , F ≥ 0 (15) or is equal to η ^vep. Advection of material is performed using a marker-and-cell technique (Gerya and Yuen, 2003) and we compute rigid body rotation of stresses analytically (see also Candioti et al., 2021, for a more detailed description of the algorithm). A marker chain and a stabilisation algorithm (Duretz et al., 2011b, 2016) is employed to allow for dynamic evolution of topography. Surface processes, namely erosion and sedimentation, are parameterized in the model. In case, the newly computed topography falls below a sedimentation level of z = −5 km it is corrected back to the sedimentation level. The resulting gap between the newly computed topographic level and the sedimentation level is filled instantaneously with markers describing a calcite or a mica rheology. The rheology of the sediment markers alternates every 2 Myr. In case the newly computed topography rises above an erosion level of z = 2 km, the topographic level is corrected downward. The distance of correction is calculated by multiplying the current time increment, Δt, by a constant erosion velocity of 0.5 mm yr⁻¹. All material parameters used in the models presented here are listed in Table 1.

Here, F D = 2 × τ ̄ II avg is used as a representative value for the horizontal driving force per unit length. Close to the lateral model boundaries, where F _D is calculated, the smallest principal stress is vertical and is close to the lithostatic pressure, and the maximal principal stress is approximately horizontal. Therefore, F _D is a representative estimate for the plate driving force (e.g., Schmalholz et al., 2014, 2019). First, the vertical integral of τ _II(x, y) is calculated at the different horizontal positions of the numerical grid as τ ̄ II x = ∫ b a τ II x , y d y , (16) where a is the topographic level at each horizontal grid point and b = −660 km is the bottom coordinate of the model domain. The values of τ ̄ II ( x ) are then averaged horizontally inside two regions of 100 km width located at the two lateral model sides. This horizontally averaged, vertically integrated stress is termed τ ̄ II avg . The reader is referred to Candioti et al. (2020) for further detail.

We further estimate the slab-pull force per unit length by calculating the effective body force of the subducting slab. This body force, F _SP, is computed as the difference between densities of all material that is 1) enclosed horizontally within −700 ≤ x ≤ 700 km of the domain, 2) colder than T ≤ 1,300°C and 3) subducted below y = −160 km and a vertical reference density profile. This reference profile is obtained from the last time step of the cooling period of model NOSLIM by averaging horizontally all material densities at each vertical grid level below y = −160 km over the entire horizontal domain. This reference density is also used to compute the density difference, Δρ, induced by a phase transition. The body force is then calculated as follows: F SP = ∫ Ω Δ ρ g d Ω , (17) where Δρ is the aforementioned density difference integrated over the area Ω covered by material matching the three criteria defined above. This is a maximum estimate for the slab-pull force, because the surrounding mantle material exerts shear stresses on the slab surface which resist subduction. Assuming shear stress magnitudes of 1–10 MPa along the bottom and top of a ca. 500 km long slab results in force magnitudes between 1 and 10 TN m⁻¹, respectively, acting against the slab pull force due to buoyancy only. Therefore, the natural, effective slab pull force magnitudes including the shear resistance from the surrounding mantle are likely lower than the ones reported here (see also Turcotte and Schubert, 2014). However, the relative evolution of F _SP likely remains as reported here.

Convergence is started at 110 Myr in model history. Subduction initiates at ca. 113 Myr at the transition zone between the proximal and the necking domain within the mantle lithosphere of the right margin (Figures 3A,E). SZI occurs at a region where no serpentinite is present which indicates that the initially horizontal serpentinite layer is not important for SZI. This conclusion is also supported by the findings of Candioti et al. (2021), who report successful SZI with and without a serpentinite layer in similar numerical models. Maximal bending stresses within the lithosphere of the subducting and overriding plate are ca. 500–700 MPa and limited to a 10–25 km thin layer (Figure 3E). Until ca. 129 Myr, large portions of the serpentinite are sheared off the subducting slab (Figure 3B). The values for stresses decrease below ca. 100 MPa within the overriding plate and remain at ca. 250–400 MPa where the subducting plate is bent below the overriding plate (Figures 3F,G). This stage marks the successful SZI in this simulation. The embryonic ocean basin is closed at ca. 140 Myr in model history. At this stage, the serpentinite material has largely been sheared off the subducting plate and reorganized along the plate interface (Figure 3C). A thrust wedge has formed at ca. 160 Myr and a large part of the subducting slab has been detached from the subducting plate (Figure 3D). The values of deviatoric stresses within the thrust wedge remain between ca. 100–250 MPa (e.g., − 200 ≤ x ≤ −100 km and y > − 30 km in Figure 3H).

Figure 4 shows the stage of SZI in the different models. In the reference model, NOSLIM, and in SLIM300 a shear zone forms below the right margin (blue region in right column of Figures 4A,B). Values for dissipative heating inside the shear zone are above 10⁻⁵ W m⁻³ (red contour line in Figures 4A,B) and the temperature inside the evolving shear zone increases compared to the surrounding ambient mantle temperature. This temperature increase is indicated by the geometry of the isotherms which are in shallower depth inside the shear zone compared to the corresponding isotherms to the left and right side of the shear zone. By comparing the vertical position of the grey isotherms in Figures 4A,B, one can estimate the temperature rise inside the shear zone assuming that before shear zone formation the isotherms were approximately horizontal, like in the corresponding left margins. For example, at the right margin of model NOSLIM, the isotherm for 727°C inside the shear zone is in approximately the same depth as the horizontal isotherm for 637°C to the left of the shear zone. Therefore, in NOSLIM the temperature rises by ca. 90°C and in SLIM300 by ca. 50°C, using the same estimate based on isotherm geometries. Kiss et al. (2019) showed that a small temperature increase of only ca. 50°C already can cause strain localization by thermal softening. Compared to models NOSLIM and SLIM300, less dissipative heat is generated by the shear zone forming in model SLIM200. The temperature rise inside the shear zone of model SLIM200 is still ca. 50°C, which is comparable to model SLIM300. In SLIM100 the heat generated by the conversion of mechanical work into heat is insignificant and the temperature inside the shear zone rises by less than 10°C (compare Figures 4A,B to Figures 4C,D). Figure 5 shows which regions in the model exhibit only rate–independent plastic, or brittle, deformation (highlighted in blue) and which regions exhibit elastic or rate–dependent plastic, or ductile, deformation (highlighted in white, see Section 2.3 for calculation). In contrast to the ductile (white regions in Figures 5A,B) shear zones formed in NOSLIM and SLIM300, the shear zones in SLIM200 and SLIM100 are mainly formed due to brittle deformation (blue regions in Figures 5C,D). The reason is that in SLIM200 and SLIM100 the stress limits are lower and, hence, the regions in which stresses exceed the stress limit (i.e., here the brittle region) become larger. Further, brittle shear zones also form below the left margin in SLIM100 (Figure 4D and Figure 5D). In consequence, the subduction zone in SLIM100 evolves completely different compared to the subduction zone in SLIM300. Subduction in SLIM100 is initiated below the opposite margin (Figures 6A,B) and the slab is subducted much slower (compare depth of slab tip in the left and right column of Figure 6C). The detachment of the slab occurs approximately at the same depth in SLIM100 and SLIM300, but much later in SLIM100 compared to SLIM300 (156 and 150 Myr, respectively, Figure 6D). We consider the slab as detached when the distance between the 1,300°C isotherms bounding the slab is less than 10 km. The slab in model SLIM300 is much longer than it is thick, thus resembling the typical rectangular cross-sectional shape of a slab. In contrast, the shape of the slab in SLIM100 resembles more the cross-sectional shape of a lithospheric droplet rather than a subducting slab (compare left and right panel Figure 6D). The deformation behaviour observed in SLIM100 is very similar to the so-called drip-off mode described in the numerical simulations of Thielmann and Kaus (2012). This deformation behaviour does not correspond to the subduction of a coherent plate-like slab.

As shown above, the maximum temperature increase associated with the formation of shear zones is less than 100°C. It is, therefore, difficult to isolate the temperature increase exclusively caused by shear heating because of the simultaneously ongoing temperature advection and conduction. To clearly show that the shear zones, causing SZI in our models, are due to shear heating and associated thermal softening, we have performed two additional simulations of NOSLIM and SLIM300 in which we have deactivated the shear heating (Figure 7). This deactivation was done by setting the term H _D in Eq. 3 to zero. When shear heating is deactivated (see right column in Figure 7) SZI is not successful. Instead, buckling and thickening in the regions of the two margins is the dominant deformation behaviour of the lithosphere. Therefore, the four simulations NOSLIM and SLIM300 with and without shear heating clearly show that shear heating is causing SZI in our models. Consequently, other model features such as the mechanical heterogeneities, mimicked by initial elliptical heterogeneities, or the serpentinite layer are not able to cause SZI.

In Figure 8 we report the evolution of the plate driving forces per unit length, F _D, and an approximate estimate for the slab pull force per unit length, F _SP, over time. The magnitude of F _D increases to ca. 18 TN m⁻¹ in NOSLIM, ca. 15 TN m⁻¹ in SLIM300 and ca. 12 TN m⁻¹ in SLIM200 between the onset of convergence and ca. 113 Myr (Figure 8A). A maximum F _D = 14 TN m⁻¹ is reached at ca. 122 Myr in SLIM200. Between ca. 113 and 114 Myr, values for F _D decrease in NOSLIM and SLIM300 and remain relatively constant in SLIM200. In NOSLIM, SLIM300 and SLIM200, the values for F _D overall decrease to ca. 10 TN m⁻¹ between ca. 130 Myr and ca. 135 Myr. In SLIM100 the magnitude of F _D continuously increases to ca. 10 TN m⁻¹ until ca. 122 Myr in model history (blue solid line in Figure 8A).

The density within the subducting slab increases by up to 150 kg m⁻³ compared to the ambient mantle density structure when it transects y = −410 km (Figures 8D–F). The reason for this increase in density is the olivine–wadsleyite phase transition. Due to a positive Clapeyron slope (i.e., the change in temperature–pressure space), this phase transition occurs at shallower depths for cold material than for hot material. Therefore, the phase transition is deflected upward inside the cold subducting slabs (see Figures 8D–F) and induces an additional pull force. A gradual increase in F _SP until the onset of necking in the subducting slab (left edge of grey rectangles in Figure 8B, see also stages in Figures 8D–F) is observed in all models. This stage is reached after ca. 141 Myr in SLIM200, ca. 150 Myr in SLIM300 and ca. 158 Myr in NOSLIM. The corresponding magnitudes are F _SP ≈ 15 TN m⁻¹ in SLIM200, F _SP ≈ 25 TN m⁻¹ in SLIM300 and F _SP ≈ 35 TN m⁻¹ in NOSLIM. When the slab detaches and sinks into the mantle, magnitudes of F _SP increase rapidly to ca. 40 TN m⁻¹ in SLIM200, ca. 42 TN m⁻¹ in SLIM300 and again ca. 40 TN m⁻¹ in NOSLIM. This rapid increase of F _SP occurs ca. 4, 13 and 21 Myr after the slab has transected y = −410 km in SLIM200, SLIM300 and NOSLIM, respectively. The depths of slab detachment also vary from y ≈ − 200 km in NOSLIM, y ≈ − 150 km in SLIM300 to y ≈ − 30 km in SLIM200 (Figure 9). In SLIM100 the magnitude of F _SP increases slowly to ca. 30 TN m⁻¹ until ca. 156 Myr. At this stage, drip-off of the lithosphere occurs in this model (Figure 6D).

In combination with temperature-dependent rheologies, shear heating is potentially a feasible mechanism to form a shear zone transecting the lithosphere (Yuen et al., 1978; Thielmann and Kaus, 2012; Jaquet and Schmalholz, 2018) and initiate subduction at a hyperextended continental margin (Candioti et al., 2021). Kiss et al. (2019) derived an equation to predict the maximum temperature rise inside a shear zone, T _SH (Eq. 22 in Kiss et al., 2019), that is formed via shear heating. This prediction relies only on material parameters and the applied absolute horizontal velocity difference, and has already been applied to SZI at homogeneous hyperextended continental margins (Kiss et al., 2020). According to their prediction, a temperature rise by 50°C is sufficient for spontaneous shear localisation via shear heating and the associated thermal softening. In our models, the dissipation created by the conversion of mechanical work into heat is ca. 10⁻⁵ W m⁻³. Consequently, the temperature increases by 50–90°C in models SLIM200, SLIM300 and NOSLIM. Due to thermal diffusion such a moderate temperature increase is likely not detectable in natural shear zones having a thickness in the order of a few kilometers (e.g., Kiss et al., 2019). The temperature inside the evolving shear zone can be predicted by the equation of Kiss et al. (2019) (deflection of grey solid lines in Figures 4A,B). Therefore, shear heating is efficient for stresses inside the lithospheric mantle that are as low as ca. 200 MPa in our models. This stress magnitude represents the maximum stress during the onset of the formation of a shear zone and stresses are decreasing during progressive shear zone formation (Kiss et al., 2019). This observation is crucial as shear heating is often critized for requiring too high stresses to be an efficient mechanism for ductile shear zone localization (e.g., Platt, 2015). Furthermore, Kiss et al. (2019) reported that shear heating is already efficient for stresses as low as 200 MPa. The predicted temperature inside the shear zone forming in SLIM200 is still accurate, although the deformation inside and around the shear zone is mainly occurring in the brittle-plastic field due to the low stress limit (Figure 5C). Although the generated dissipation is relatively lower compared to higher stress limiters, the temperature increase in SLIM200 and SLIM300 are comparable (compare isotherms in Figures 4B,C). The estimated temperature increase of ca. 50°C in SLIM200 is the limit for shear heating causing spontaneous shear localisation predicted by Kiss et al. (2019). This indicates that likely ductile and brittle-plastic processes are equally important for shear zone localization in model SLIM200. In case stresses are lower than 200 MPa, brittle-plastic processes dominate the deformation inside the lithospheric mantle and shear heating is likely inefficient (flat grey isotherm in Figure 5D). Thielmann and Kaus (2012), Jaquet and Schmalholz (2018) and Kiss et al. (2020) already showed that in absence of other softening mechanisms, thermal softening is required for convergent plate boundary formation and successful SZI in lithospheric-scale geodynamic models. We here confirm these results. In the models we present, thermal softening is indeed the responsible mechanism that causes SZI (see Figure 7).

We report a minimum magnitude of ca. 14–15 TN m⁻¹ for horizontal plate driving forces required for SZI via shear heating (SLIM200 and SLIM300, Figure 8A). While a magnitude of ca. 30–37 TN m⁻¹ for plate driving forces was required to initiate subduction at a homogeneous continental margin in the models of Kiss et al. (2020) and Auzemery et al. (2021), we here report a maximum value of ca. 18 TN m⁻¹ for a continental margin that includes elliptical mechanical heterogeneities. Subduction initiation at ca. 10 TN m⁻¹ is observed in SLIM100. However, subduction initiation in this model is controlled by brittle-plastic deformation rather than by shear heating. In our model SLIM100, the lithospheric deformation style changes from a subduction (Figure 6C) into a delamination, or drip-off, style with ongoing convergence (Figure 6D). Such drip-off deformation behaviour is consistent with earlier studies on subduction initiation within mechanically weak lithospheric plates (e.g., Thielmann and Kaus, 2012). We therefore propose that the minimum driving force magnitude predicted by SLIM200 and SLIM300 of 14–15 TN m⁻¹ are likely more meaningful estimates for SZI at natural geodynamic settings than the magnitude predicted by SLIM100. This magnitude could likely be further reduced by modulating, for example, the olivine dislocation and Peierls flow law parameters.

The density difference between the slab and the surrounding mantle due to the uplift of the olivine–wadsleyite phase transition (at y = −410 km) in our models is between ca. 150–200 kg m⁻³. This petrological-thermomechanically coupled density prediction is slightly lower but overall in good agreement with previous estimates (275 kg m⁻³, Turcotte and Schubert, 2014). This increase in density causes an increase in slab pull force magnitude of ca. 16 TN m⁻¹ (Turcotte and Schubert, 2014). Although this estimate is based on buoyancy only and does not consider shear forces along the subducting slab, it is meaningful as these shear forces already act on the slab before the phase change. Our model SLIM200 confirms this estimate. A sharp increase in slab pull force occurs when a magnitude of ca. 17 TN m⁻¹ is reached at ca. 141 Myr (Figure 8B). In SLIM200, SD begins immediately after the slab has transited a depth of 410 km (Figures 8C,F). This observation indicates that slabs must at least be strong enough to support the additional pull force induced by the phase transition in order to penetrate further into the mantle without detaching immediately after the olivine–wadsleyite phase transition occurs within the slab.

In our models we observe the following correlation between the slab strength and the timing of slab detachment: The weaker the slab, the earlier it detaches after transiting the olivine–wadsleyite phase transition. We further observe a correlation between the strength of the slab and the depth of the detachment: weak slabs detach already in the horizontal region of the subducting plate, whereas strong slabs detach in ca. 100–200 km depth (Figure 10). In all models, the slab detaches within ca. 1 Myr from the onset of necking. These observations are consistent with earlier studies of slab detachment (e.g., Andrews and Billen, 2009; Baumann et al., 2010; Duretz et al., 2011a; Schmalholz, 2011; van Hunen and Allen, 2011; Thielmann and Schmalholz, 2020). The estimates for slab pull forces we report here are between ca. 18 and 35 TN m⁻¹. These magnitudes are lower than previous calculations (e.g., Forsyth and Uyeda, 1975). However, these magnitudes only consider forces arising from buoyancy contrasts between the subducting slab and the surrounding mantle and do not include viscous tractions along the slab. Such tractions reduce the effective slab pull force (e.g., Schellart, 2004).

Natural orogens like the Variscides, the Pyrenees or the European Alps may have resulted from the closure of narrow oceanic basins (e.g., Chenin et al., 2019). These immature extensional systems are characterized by insignificant production of new, mature oceanic lithosphere and subduction of these embryonic oceans produce insignificant amounts of volcanic arc–basalts (e.g., McCarthy et al., 2020). In absence of a significant active mid–ocean ridge, such embryonic basin–margin system is likely more laterally homogeneous (due to lack of a ridge) compared to an ocean basin including a significant spreading ridge. Consequently, embryonic ocean basins likely lack significant lateral variation in temperature and, therefore, strength. Hence, the closure of such basins from SZI to SD may be considered as, respectively, natural compressive and tensile failure tests of the lithosphere. The models we present have significant implications for these natural systems. For example, subduction in the Western Alpine-Tethys presumably initiated below the hyperextended continental margin of the Adriatic plate at ca. 85 Ma (e.g., Manzotti et al., 2014) and lead to closure of the embryonic Piemonte–Liguria ocean (e.g., McCarthy et al., 2020). Our models predict that, compared to a weak subducting lithosphere, a strong lithosphere requires higher forces to fail under compression and yield SZI and plate boundary formation. However, under tension, which develops due to increasing pull forces induced by the subducting slab, the same strong lithosphere may be able to hold a continuous subducting slab down to 660 km depth over more than 20 Myr after basin closure. In contrast, a weak lithosphere fails immediately under tension and exhibits shallow SD early after closure of the embryonic ocean basin. This result may help explaining the strong east–west variability of slab geometries observed in tomographic images from the present-day European Alps (e.g., Lippitsch et al., 2003; Hua et al., 2017; Kästle et al., 2019).

Our models NOSLIM and SLIM300 (Figures 9A,B) suggest that the embryonic Piemonte–Liguria ocean might have been weaker in the east compared to the centre which lead to early and shallow slab break-off in the eastern Alps (slab gap, Lippitsch et al., 2003) whereas the slab remained more or less intact in the central region (Zhao et al., 2016). The resulting asthenospheric window in the slab (see glyphs in Figure 9B) below the eastern Alps might have induced a heat pulse in the crustal units (e.g., Davies and von Blanckenburg, 1995). For example, in simulation SLIM300 the 1,300°C isotherm is uplifted towards shallower depth in the region of SD (Figure 9B). The 1,300°C isotherm is uplifted to approximately the same depth as that of the horizontal 800°C isotherm to the left and right side of the region with SD. Therefore, the temperature can be locally increased by ca. 500°C with respect to temperatures in an undeformed lithosphere. Such local temperature increase due to SD might contribute to magmatism associated with the Bergell intrusion, which has been explained by SD (e.g., Davies and von Blanckenburg, 1995). Certainly, the Alps have a distinct 3D architecture and our models are only two-dimensional. Previous three-dimensional modelling studies have studied the impact of along-strike variations on the process of slab detachment (van Hunen and Allen, 2011; von Tscharner et al., 2014; Duretz et al., 2014). Because of their two-dimensional nature, our models do not capture along-strike propagation of slab detachment. Nevertheless, they successfully account for the physics of necking which is the underlying physical process governing slab detachment (e.g., Schmalholz, 2011; Duretz et al., 2012). Moreover, lateral propagation of slab detachment is triggered by asymmetries in the lateral structure of collision zones (van Hunen and Allen, 2011), which remain difficult to constrain through the course of the Alpine history. We hence suggest that our models remain applicable, at first order, to the study of slab detachment in the Alps.

Slab detachment is a short-lived process and we just mention also briefly two regions in which SD is presumably occurring at present day. In the Vrancea area (e.g., Wenzel et al., 1998) a nearly vertically hanging slab which is still attached to the overlying plate induces earthquakes at an intermediate depth between 70 and 180 km. However, the maximum depth of the slab might be much deeper (Sperner et al., 2001; Ferrand and Manea, 2021). Also, recent seismic data from the Hindu Kush indicates a potential onset of slab detachment or delamination at depths of ca. 200 km (e.g., Kufner et al., 2021). Such detachment depths are well predicted by our models (Figures 8D–F).

Three-dimensional numerical models of global mantle convection have been used to constrain the effective strength of the lithosphere. Mallard et al. (2016), for example, applied a constant von-Mises stress limiter function in their models (similar to this study) and showed that stress magnitudes between 150 and 200 MPa in the lithosphere result in the most realistic plate-like behaviour. The model results of SLIM200 confirm the results reported by Mallard et al. (2016). However, our models predict plate-like behaviour also for stresses higher than 200 MPa in the lithosphere (see for example model NOSLIM), likely because of the higher numerical resolution. In our models, high stresses in the lithosphere are limited to a 10–25 km thin, load-bearing layer (see f.e. Figures 3E–H). Such thin layers cannot be resolved numerically in most global mantle convection models. In comparison, dividing, for example, F _D = 18 TN m⁻¹, which was necessary to overcome the plate strength in NOSLIM and initiate subduction, by a ca. 50 km thick lithosphere yields an average deviatoric stress of ca. 180 MPa which agrees with the equivalent von-Mises shear stress between 150 and 200 MPa of Mallard et al. (2016), which is required to generate a plate size–frequency distribution observed for Earth. Therefore, the magnitude of the horizontal driving forces in our numerical simulations agree with horizontal forces predicted by recent global mantle convection models.

Because of the high numerical resolution, our models can test first the resistance of plates to failure under variable applied compressive stresses (during SZI) and subsequently under variable tensile stresses (during SD). Therefore, models as presented here may represent a complementary method to classical laboratory experiments, inversion from geophysical data and global mantle convection models and help to better constrain the effective strength of the lithosphere.