Influence of magma-poor versus magma-rich passive margins on subduction initiation

to early that, the process, may be discovered which could affect the and all legal disclaimers that to the pertain. the oceanic lithosphere. Abstract We present a new numerical modelling study of subduction initiation at (hyper-extended) magma-poor and magma-rich continental passive margins. In particular, we test how the structure and rheological stratification of these two end-member types control the formation and thermo-mechanical evolution of subduction zones. The serpentinization of mantle lithosphere in a magma-poor continental rifted margin leads to rheological decoupling at the base of the continental crust, which induces shear localization during subsequent shortening. Under these conditions, a shear zone propagates into the mantle lithosphere and leads to subduction initiation at the transition between ocean and passive margin. In contrast, a magma-rich rifted rheologically coupled, which creates a buttressing effect and transfer of deformation into the oceanic domain during shortening, where the subduction zone initiates. These results are quantitatively compared and in agreement with the geological record of subduction initiation in the Alps and Dinarides – Hellenides, where the relics of end-member types of continental rifted margins of the same Adriatic micro-continent bordering the Alpine Tethys and Neotethys oceans, respectively, are observed.


Introduction
Recent studies have shown that the mechanics of subduction initiation remains a critical plate tectonic process that is less understood (Hall, 2019;Stern & Gerya, 2018). These studies emphasize that pre-existing tectonic plates configurations and density contrasts across mechanically weak zones are key ingredients for subduction initiation (e.g. Arcay et al., 2020;Nikolaeva et al., 2010). These conditions can be met at oceanic transform and major fracture zones or in mantle plume settings, where strong lateral density and strength contrasts exist (e.g. Cloetingh et al., 2021;. Oceancontinent transitions also meet these conditions, but direct evidence for subduction initiation at passive margin is scarce (Crameri et al., 2020;Gurnis et al., 2004). These observations and models have led to the idea that externally driven processes such as subduction zone invasion are needed for the genesis of subduction zones at continental margins (Crameri et al., 2020;Duarte et al., 2013;Zhou et al., 2020). Although this mechanism may explain some natural examples (e.g. western edge of the Pacific plate, Hall, 2019;South Sandwich, van de Lagemaat et al., 2021), its spontaneous retreating mode of subduction (Zhou et al., 2020) contrasts with the geological record of subduction initiation at magma-poor continental margins such as in the Pyrenees or the European Alps orogenic systems (McCarthy et al., 2018), suggesting that subduction initiation must have occurred at passive margins in the plate tectonics history.
These modelling studies have shown that, under compressive stress regime, deformation is focussed in the ductile lower continental crust when the oceanic lithosphere at the margin is older than ~50 Myr and the continental lithosphere is relatively thin (~100-150 km). When such conditions are met, subduction initiation is physically feasible by the evolution of major shear zones in the ductile crust that propagates into the mantle where thermal softening facilitates the formation of a subduction zone (Auzemery et al., 2020;Kiss et al., 2020). In contrast, deformation localises within the rheologically weaker oceanic lithosphere, when they are younger than 50 Myr and/or in case of high convergence rate (> 2 cm yr -1 , Zhong & Li, 2019) creating intra-oceanic subduction.
All previous modelling studies have not yet quantified the impact of significant differences in geometry and rheology between magma-rich and magma poor continental passive margins (MRPM and MPPR, respectively) for subduction initiation (Fig. 1), constrained by the available observational record (Geoffroy, 2005;Sengör & Burke, 1978;Tugend et al., 2020;White et al., 1987).
In particular, MRPM are typically driven by mantle upwelling and are associated with significant magmatic activity during rifting (Fig. 1a, Corti, 2009;Li et al., 2020;White et al., 1987). Important is the presence of magmatic intrusions and extrusions as well as underplating of mafic bodies at the base of the crust, which significantly increases its effective strength after cooling (Kelemen & Holbrook, 1995;Paton et al., 2017). In contrast, MPPM are typically driven by lithospheric stretching without significant decompressional melting magmatism, resulting in exhumation and hydration of the mantle in the footwall of major detachments (Fig. 1b, Brun & Beslier, 1996;Clerc et al., 2018;Manatschal et al., 2001;Mohn et al., 2012). The hydration-induced serpentinization penetrates into the mantle lithosphere and creates rheologically weak decollement layers that are prone to tectonic reactivation (Guillot et al., 2015;Hilairet et al., 2007).  Geoffroy, 2005); and b) magma-poor continental rifted margins (inspired from Mohn et al., 2012) for two evolutionary stages.
In this study, we aim to understand and quantify subduction initiation in the vicinity of contrasting end-member types of passive continental margins (MRPM and MPPR) during far-field induced shortening by means of thermomechanical numerical modelling, constrained by the available observational record. In particular, we test the influence of differences in crustal and lithospheric structure, composition and evolution of deformation. A two-steps strategy is adopted to understand both first order controlling mechanisms and the more detailed response. In a first step, we assess the deformational response to the degree of rheological crust-mantle coupling at continental passive margins. In a second step, we use more realistic geometries and compositions of magma-poor and magma-rich continental passive margins to test their specific impact on subduction initiation. The modelling results are validated against observations and quantitative reconstructions of the Adriatic micro-continent, an ideal well-studied area that contains the record of both endmember types of passive continental margins (MRPM and MPPM). This record is exposed in two European orogenic systems, the Dinarides-Hellenides and the Alps, for the MRPM and the MPPM, respectively.

Methods
We performed two-dimensional (2D) thermo-mechanical modelling to simulate the shortening of a visco-elastic-plastic lithosphere. The finitedifference, marker-in-cell code (MDoodz; Duretz et al., 2016) solves the equations of momentum (1) and mass conservation (2), as well as the heat equation (3) formulated as: where v is the velocity vector, T is the temperature, k is the thermal conductivity, is the density, C p is the heat capacity, Q r is the radiogenic heat production, is the deviatoric stress tensor, is the deviatoric strain rate tensor, P is the pressure and g is the gravity acceleration vector. The term SH = vp describes the production of heat by visco-plastic (vp) dissipation (shear heating).
For details regarding the mathematical model and algorithms, see supplementary material S1.
The model domain is a section of 3000 × 500 km and the grid resolution is 1 x 1 km in both dimensions (Fig. 2). The model top boundary is a true free surface (Duretz et al., 2016), computed as: where is the height of the free surface, is the topographic diffusivity and ℎ is the sedimentation rate. Because the effects of sediments have been inferred to be critical for the development of a self-sustaining subduction zone (Cloetingh et al., 1984;Nikolaeva et al., 2010;Sobolev & Brown, 2019), erosion and sedimentation are parameterized in the equation, 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. To limit the loading effect, we adopted a low sedimentation/erosion rate, setting diffusivity to 10 -6 m 2 .s −1 (for h > 0), and sedimentation to 0.5 mm.yr −1 (for h < 0). In all models, the accommodation space is filled with a sedimentary material composed of alternating calcite and mica layers, which is a rheologically weak composition.
A constant normal velocity is applied at the left and right boundary of the model (|+V in | = 0.625 cm.yr -1 ) from the surface to a depth of 200 km. The total convergence rate of 1.25 cm.yr -1 reflects an average rate of slowly converging systems, in accordance with estimates specific to the closure of the Vardar/Neotethys and Piemont Liguria Ocean (e.g., Capitanio & Goes, 2006;Handy et al., 2010;Le Breton et al., 2020). In order to satisfy mass conservation, an outflow velocity (-V out ) is distributed at the base of the model and on the sides White and black arrows show the directions of outflow (-Vout) and the velocity applied at the left boundary of the models (+Vin), respectively. PM is the passive margin.
from 200 km depth to the base of the box (Fig. 2). The outflow is proportionally distributed over the length of each boundary.
Deformation is governed by Drucker-Prager plasticity, dislocation creep, diffusion creep, Peierls creep and elasticity (see parameters in Tables 1 and 2 and supplementary material S1 for equations). Following previous studies Kiss et al., 2020;Zhong & Li, 2019) we account thermal softening instead of pre-defined strain softening mechanisms to allow for shear localization leading to subduction initiation. Such a setup avoids boundary conditions artefacts during modelling. We however acknowledge that fluid driven softening mechanism may also be critical for intra-oceanic subduction initiation (Dymkova & Gerya, 2013).
Ref.  Here h is the thickness of the different layers for the ocean (O) and the continent (C), ρ is the density, k is the thermal conductivity, Qr is the radiogenic heat production, ϕ is the friction angle, A is a pre-factor, V is the activation volume, f is a correction factor, n is the stress exponent, G is the shear modulus, Q is the activation energy. The cohesion C is set to 20 MPa for all lithologies. References for rheology are R_95: Ranalli (1995); H&C_83: Hansen & Carter, 1983;H&K_03: Hirth and Kohlstedt (2003); M_98: Mackwell et al., 1998;R&D: Rybacki & Dresen, 2004;K_90: Kronenberg et al. (1990); H_07: Hilairet et al. (2007).

Modelling strategy
We have chosen a two-steps modelling strategy: (1) we initially run a series of reference models consisting of a simplified rheological stratification where the entire crust is represented by the same material, and (2), we perform another set of models that assumed a layered and laterally variable crustal composition that matches as close as possible the two end-member types of passive continental margins (MPPM & MRPM).
The models (Fig. 2) consist of a continental lithosphere made of a 36 km thick crust and a lithospheric mantle with the thermal base (H L ) is fixed at a depth of 170 km reflecting a lithosphere of Palaeozoic thermo-tectonic age (Artemieva, 2009;Burov, 2011 (Parsons & McKenzie, 1978;Houseman et al., 1981;Turcotte & Schubert, 2014). Such a configuration, together with a relatively slow convergence rate favours subduction initiation at passive margins (Auzemery et al., 2020;Kiss et al., 2020;Zhong & Li, 2019). The passive margin geometry includes a crust that wedges-out towards the ocean-continent transition. The base of the lithosphere (controlled by the 1330°C isotherm) at the passive margin is defined by a linear interpolation between the ocean and the continent.
The initial geotherm is computed assuming steady state and accounts for radiogenic heat production in each layer and a constant asthenosphere temperature of 1330°C. The resulting temperature at the continental Moho is ~460°C. When displaying modelling results, we focus on the mechanically significant part of the lithosphere defined by the 1000°C geotherm (~base of the strength envelope, Burov, 2011).
In this first step, we carried out 32 numerical models and systematically varied the length and the crustal composition of the passive margin lithosphere, thereby reducing the mechanical coupling between the crust and the mantle lithosphere at the margin. The reference models test 8 different lengths of passive margins changing from 50 km to 400 km, covering most of natural passive margin lengths (Reston & Manatschal, 2011), as well as 4 different rheologies for the whole crust (Table 1). The modelling results are used for constructing domain diagrams (Figs. 3a, b), which allow for mapping the mode of deformation against the length and the crust-mantle coupling of the margin.
They thus provide an overview of suitable parameter combinations for subduction initiation and provide a geometric and rheological frame for models with specificities tailored to represent magma poor and magma rich passive margins.
The second set of models includes the specific crustal characteristics of MPPM and MRPM as illustrated in figure 1, with the same thermal base in the oceanic and continental domain for both models. A different flow law was attributed to each material, in agreement with the simulated type of oceanic or continental layered lithosphere, magmatic material and overlying sediments, as well as their lateral distribution, to account for a representative simplification of a more realistic natural scenario (see Figs. 6a, 7a). To be able to highlight the relevance of crust-mantle coupling and to correlate models with our parametric study, the crust and mantle thicknesses were kept similar. However, a 2-layers model is now used for the continental crust consisting of an upper crust of 20 km thickness with dry quartzite rheology and a lower crust between 20 and 36 km depth where deformation is described by a felsic granulite flow law (see parameters in Table 2).

Model limitation
The aim of this study is to provide insights into favourable rheological propagation of an existing subduction zone (e.g. Zhou et al. 2020). Subduction initiation in our study is driven by horizontal forcing imposed by convergent kinematic boundary conditions. Therefore, we do not account for mantle upwelling (Gerya et al., 2015;Koptev et al., 2019;Koptev et al., 2021) or mantle suction flow (Baes & Sobolev, 2017;Yamato et al., 2013). Furthermore, we do not account for climate driven variability of sediment supply to the trench controlling the lubrication of the evolving subduction zone (Sobolev & Brown, 2019). Despite these simplifications, our 2d models provide a significantly advanced approach to accurately describe the system behaviour and obtain novel insights in the key parameters controlling the mechanism associated with subduction initiation.

Reference model
Using numerical experiments, we have systematically investigated the role of the passive margin length and the degree of crust-mantle coupling on subduction initiation (Table 3 and Fig. 3a). For each rheological configuration presented above, the length of the margin (the wedge in the models) has been varied yielding three end-modes of subduction initiation: subduction at the passive margin, intra-continental subduction and intra-oceanic subduction (Figs. 4a, b, c, respectively). More specifically, intra-continental subduction refers to the subduction of the continental lithospheric mantle under the oceanic lithosphere.  Figure 3: Results of the parametric study for the reference model: (a) Area diagram showing the influence of continental crust-mantle coupling and length of the lithosphere at passive margin on the mode of subduction initiation. wgr., an., gr., dol., stand for westerly granite, anorthite, granulite, dolerite, respectively. The numbers and colours represent the experiment numbers and the modes of subduction initiation, respectively. The thickness of ductile layer at the onset of deformation is measured between the depth of the stress peak and the base of the continental crust. (b) time period of subduction initiation measured from t0 to the incipience of subduction.   (c) intra-oceanic subduction. Enlargements display the strain rate and the viscosity at the early stage of underthrusting as well as the direction of the maximum principal stress ( 1 ).

Role of crust-mantle coupling
For a constant length of the passive margin (e.g., 150 km, Fig.5), variation of continental crust rheology impacts on the degree of crust-mantle coupling and consequently controls the locus of subduction initiation. For instance, the presence of a granitic continental crust allows the decoupling of the crust from the mantle and favours subduction initiation at the margin. In this case, deformation starts at the base of the ductile crust at the ocean-continent transition, where a shear zone develops, which eventually propagates into the mantle lithosphere (Fig. 5a). In contrast, shortening of a strong doleritic or granulitic continental crust, which results in strong crust-mantle coupling, consistently leads to folding of the lithosphere followed by intra-oceanic subduction (Figs. 5c,d). In this case, early-stage deformation affects the passive margin, which bends downward (Fig. 5c). However, strain does not localize over a long period of time and deformation shifts from the passive margin towards the oceanic domain where subduction occurs.
Models where continental crust is simulated using an anorthite flow law, (intermediate case of coupling) follow a similar evolution as those with a granitic crust. However, subduction at the margin is intricate and folding is observed in the oceanic domain (Fig. 5b). We note that folding is better expressed close to the passive margin than in the middle of the oceanic plate (Figs. 5b, c), likely reflecting stress concentration at the ocean-continent transition. Evidence of this mechanism is highlighted by a decrease of the amplitude and the wavelength of folding from the margin to the ocean (Fig. 5b). In case of stronger coupling (dolerite crust), deformation at the margin is limited and folding mainly affects the oceanic domain (Fig. 5d). Such early intra-oceanic folding leads to higher sedimentation in the newly formed flexural basins (Fig. 5d).

Role of the passive margin length
Our results show that subduction initiation at passive margins requires the development of a shear zone that subsequently propagates into the mantle lithosphere (Fig. 4a). A small wedge angle represents a case of a long passive margin where the ductile crust gets gradually thinner towards the ocean-continent transition.
When the thickness of the ductile crust gets lower than ca 3.5 km, stresses within this layer get too high and the ductile crust does not act as a decollement anymore. Instead, deformation is transmitted further into the ductile crustal layer toward the continent (Fig. 4b). In this situation, shortening is accommodated by folding of the mantle lithosphere at the margin, resulting in deflection of the continental Moho. In the later stage of the model, strain propagates into the mantle leading to intra-continental subduction (Fig. 4b). The oceanward polarity of the subduction zone is probably controlled by the buttressing effect of the long and overall strong passive margin crust, which is coupled to the mantle lithosphere in the area of the ocean-continent transition. Therefore, our models suggest that the formation of a subduction zone at passive margins is more easily achieved with a steep and decoupled passive margin (Fig. 3a). In other cases, deformation is transmitted toward the continental crust (Fig. 4b) or the oceanic lithosphere (Fig. 4c).

Timing of subduction initiation
Our results highlight an important correlation between geometrical (length of passive margin and margin angle) and rheological (thickness of ductile continental crust) initial conditions with the period of time required for subduction initiation (Fig. 3b). We refer to duration of subduction initiation in our models as the period between the start of shortening (at t=0) and the moment when the continuous shear zone analogous to an incipient subduction plate boundary is established. Furthermore, we use the thickness of the ductile crust as a proxy for the degree of crust-mantle coupling, where a thicker ductile crust represents a decoupled case and vice versa. For a similar convergence rate, initiation of subduction requires less time for a decoupled lithosphere (<10 Myr) than for a coupled lithosphere (>15 Myr, see Fig. 3b), because a weak lower crust favours strain localization early in the evolution of a subduction zone. For the same reason, the length of the margin increases the duration of subduction initiation.
An exception exists with strong crust mantle-coupling at the passive continental margin (e.g., crust simulated using dolerite). In such cases subduction initiation does not occur at the passive margin, which deforms over a long period of time, but within the oceanic plate (Fig. 5d). We emphasise that the duration of subduction initiation also includes early phases of distributed deformation of the passive margin prior to strain localization. This phase is longer by ca. 10 Ma for weak margins (e.g., anorthite crust) as opposed to strong margins consisting of doleritic crust as indicated in figure 3b.

Application to magma-poor and magma-rich passive margins
The parametric study allowed us to infer the sensitivity of parameters controlling localization of deformation and subsequent subduction initiation. The next step is to understand the impact of a more realistic 150 km long magmapoor ( Fig. 6) and magma-rich (Fig. 7) passive margin geometries and composition for subduction initiation.

Subduction initiation at magma-poor passive margins
The magma-poor passive margin model (Fig. 6a) accounts for the necking of continental crust, resulting in a large omission of the lower crust in distal margin domains and the exhumation and hydration of sub-crustal mantle peridotites at the ocean-continent transition (Brun & Beslier, 1996;Manatschal et al., 2001;Mohn et al., 2012). As in geophysical models, the exhumed mantle shows variable degrees of serpentinization (Bayrakci et al., 2016), which is accounted for in the model by placing distributed patches of serpentinite material within the uppermost mantle lithosphere, where the mantle is located at the sea floor and below the extended margin (Fig. 6a). The latter represents serpentinization due to normal faulting cutting across the extended margin crust.
By implementing these complexities, we account for a relatively weak lithosphere at the ocean-continent transition that represents a thinned continental lithosphere overlying a serpentinized mantle, and therefore, a relatively weak crust-mantle coupling (Fig. 6a).
The evolution of our simulation follows the general trend of a model with crust-mantle decoupling at the base of the margin (Fig. 6b)   The results show that the thrust-type deformation sequence is younging towards the continent (fault 1 is older than fault 2, Fig. 6b) during the early stage of deformation. Different to the reference models, strain rates are highest (~10 -13 s -1 ) within the serpentinized mantle, (Fig. 6b1) and not within the ductile crust.
The model shows that at ca. 10 Ma of shortening, shear along the weak serpentinites links up with shear within the ductile crust below the necking zone to form a throughgoing shear zone at the base of crust (Fig. 6b2). Furthermore, this shear zone is linked to a major thrust-type structure (number 2 in Fig. 6b2) that evolves near the necking zone where the margin is thickening, and the wedge angle increases. Subsequently, the distal margin is thrust below the continental lithosphere and accreted to the overriding plate at depth of 25-35 km ( Fig. 6b3). At the position of the former distal margin a thick sedimentary wedge developed ~15 Myr after the onset of convergence.

Subduction initiation at magma-rich passive margin
The model involving a magma-rich passive margin accounts for magmatic underplating, as well as mafic intrusions and extrusions, as a result of partial melting of the mantle lithosphere (Geoffroy, 2005). Magmatic underplates are often characterized by high-velocity seismic zones in the lower crust (Kelemen & Holbrook, 1995;Sapin et al., 2021). Mafic intrusions are gabbroic complexes (sills, dikes) distributed across the MRPM. Their equivalents at the surface correspond to large basaltic lava flows identified as seaward dipping reflectors on seismic reflection data (Paton et al., 2017;White et al., 1987;). We ascribe a dolerite rheology to the magmatic underplatings while magmatic extrusions are simulated by wet olivine. In the model we implement these features by accounting for a margin that exhibits a continental crust intruded by mafic rocks.
The evolution of our model during subduction initiation of a MRPM follows the general trend of a model with a strong crust-mantle coupling at the margin (Fig.   7b). It features first distributed deformation at the ocean-continent transition, migration of deformation and folding of the oceanic lithosphere, followed by intra-oceanic subduction. The modelling shows that the addition of magmatic lithologies strengthens the passive margin close to the ocean-continent transition and increases the crust-mantle coupling, preventing strain localization. Therefore, most of the strain is localized in the oceanic lithosphere near the ocean-continent transition, while less deformation is observed in the ductile layer of the proximal continental margin (Fig. 7b1). The oceanward underthrusting takes place 350 km away from the margin, at one of the conjugate shear structures (Fig. 7b2). Although the oceanic domain is 1500 km long, strain localizes relatively close to the oceancontinent transition, similarly to the reference model (Fig. 5c). The subsequent development of a subduction zone is associated with the development of an accretionary wedge, which is recognisable ~12 Myr after the onset of convergence (Fig.7b3).

Role of the passive margin length for magma-poor and magmarich passive margin settings
We investigated the influence of margin length and associated taper angle and rheology for endmember situations of short (50 km) and long (400 km) magma rich and magma poor passive margins (Fig. 8). Overall, these models of magma-poor and magma-rich passive margins are consistent with results derived from the reference models with a granulitic and doleritic crust, respectively (Fig. 5). In the situation of a magma-poor margin, a decrease of a narrow passive margin favours subduction at the margin (Fig. 8a) similarly to the model with average (150 km) margin length (Fig. 6). In contrast, an increase in margin length does not allow strain localization at the passive margin, but at the OCT, where serpentinites act as a weak zone (Fig. 8b). In the situation of a magma-rich margin, the shortening always creates intra-oceanic subduction initiation, regardless of margin length (Fig. 8 c-d), because the mafic intrusions and underplating lead to crust-mantle coupling as well as the strengthening of the margin (Fig. 7).

Comparison to previous models
Although subduction zones often initiate through the lateral propagation from pre-existing subduction or polarity reversals (Crameri et al., 2020), other mechanisms that allow for the formation of new subduction plate boundaries must exist Crameri et al., 2020;Hall, 2019;McCarthy et al., 2018), most likely driven by far-field tectonic forcing (Gurnis et al., 2004), mantle suction flow (Baes & Sobolev, 2017), or a plume activity . In that situation, the probability of subduction initiation increases in the vicinity of continental margins (e.g. Ulvrova et al., 2019). Previous modelling studies have shown that the ability of a passive margin to convert into a subduction zone is largely governed by the age and compositional differences in the mantle lithosphere across the ocean-continent transition (e.g., Auzemery et al., 2020;Baes & Sobolev, 2017;Faccenna et al., 1999;Gurnis et al., 2004;Nikolaeva et al., 2011;Zhong & Li, 2019). Our models are in agreement with these general inferences, but further demonstrate that the degree of crustmantle coupling has also a critical impact on the localisation of subduction zones (Fig. 3), which is further controlled by the geometry and the rheological layering  (Fig. 7). Therefore, we agree that intra-oceanic subduction is more favourable in young oceanic basins (Auzemery et al., 2020;Faccenna et al., 1999;Zhong & Li, 2019), but it can also happen in old oceanic lithospheres for cases where the strength of the passive margin lithosphere is higher (e.g., through magmatism) compared to the oceanic lithosphere.
Together with the thickness of the ductile lower crust, the margin length is an important factor that influences strain localization through the crust and the mantle lithosphere. Therefore, for long (more than 300 km), decoupled and low angle passive margin, deformation is transmitted toward the continental crust ( Fig. 4b) or the oceanic lithosphere (Fig. 4c). In the latter scenario, subduction is expected to occur at a strength minima which coincides with the presence of serpentinite (Fig. 8b). Even if our concept of continental underthrusting is rarely mentioned in studies discussing plate tectonics, similar models have been obtained through reconstructions of the north African margin (Roure et al., 2012) and Cantabrian margin (Teixell et al., 2018;Viejo & Gallastegui, 2005) or in numerical studies (Kiss et al., 2020;Hamai et al., 2018).

Application to observations and reconstructions of subduction initiation near the former Adriatic passive continental margins
The evolution of continental passive margins along the former Adriatic micro-continent provides a unique opportunity to study subduction initiation in the vicinity of two continental rifted margin end-members. A MPPM and a MRPM developed over a similar pre-rifting continental structure and composition, with a similar oceanic lithospheric age at the time of subduction initiation and with comparable convergence rates. The complexity of the Mediterranean system in terms of evolution and definition of units has resulted in numerous ideas of tectonic evolution and reconstructions (Capitanio & Goes, 2006;Csontos & Vörös, 2004;Handy et al., 2010;van Hinsbergen et al., 2020;Le Breton et al., 2020;Schmid et al., 2020;Stampfli & Hochard, 2009). Our analysis of tectonic evolution ( Fig. 9) Fig.10.

Subduction initiation along the former Adriatic MPPM presently exposed in the Alps
The onset of convergence in the Neotethys Ocean was roughly coeval with the onset of Middle Jurassic opening of the Piedmont-Liguria branch of the Alpine Tethys Ocean at ~170-160 Ma (Figs. 9a-c, Liati et al., 2005;Schaltegger et al., 2015). Rifting magmatism is minor or non-existent in the Alpine chain, with the exception of the eastern Southern Alps transitional area, where a Triassic phase of extension is associated with magmatism (e.g. Lustrino et al., 2019). A MPPM has been defined across several Alpine transects (Manatschal & Bernoulli, 1999;Manatschal & Müntener, 2009;Mohn et al., 2012;Müntener et al., 2004). After 60-70 Myr of slow oceanic seafloor spreading (McCarthy et al., 2020), the closure of the Piemont-Liguria branch of the Alpine Tethys ocean was initiated by subduction at the Adriatic passive margin (Fig. 10a, Manzotti et al., 2014;Marroni et al., 2017) and was driven by the northward motion of the Adriatic plate at a convergence rate of 0.9 to 1.5 cm.yr -1 (Handy et al., 2010;van Hinsbergen et al., 2020). The Alpine Tethys Ocean was closed by the Eocene onset of continental collision, followed by the subsequent indentation of the Adriatic micro-continent in the Alps (e.g., Schmid et al., 2004a).

Figure 10:
Reconstruction of subduction initiation for the Alps and the Dinarides -Hellenides (adapted and simplified from Manzotti et al., 2014;Schmid et al., 2004bSchmid et al., , 2008. (a) Shortening initiate subduction at the Alpine Tethys Adriatic MPPM. Development of shear zones in the serpentinized exhumed mantle leads to nappe emplacement at the Adria continental margin. This evolution can be compared with our model of MPPM (Fig. 6b). (b) Development of intra-oceanic subduction in the vicinity of the Neotethys Adriatic MPPM. Locations AB, CD, of the cross-section are indicated in figures 9a and 9b. Black arrows indicate direction of extension and shortening respectively. CAB, MOR, SZZ stand for calc-alkaline basalt, mid-oceanic ridge and supra-subduction zone (magmatism).
The Adriatic passive continental margin preserved a necking domain that recorded the formation of detachments in the upper-middle crust (Mohn et al., 2012;Manzotti et al., 2014). These detachments accommodated the extreme thinning of the Adria margin and created the large tectonic omission observed between the lower continental crust and exhumation of subcontinental mantle (Manatschal & Müntener, 2009;Müntener et al., 2004). The margin had a reconstructed width of 150-250 km at the time of subduction initiation (e.g. Beltrando et al., 2010), which was controlled by the presence of shear zones in the serpentinized exhumed mantle (Beltrando et al., 2010;Malatesta et al., 2011). The subsequent nappe emplacement was associated with high pressurelow temperature (HP-LT) metamorphism (Manzotti et al., 2014 and references therein) and the development of an accretionary wedge consisting of oceanic, continental, and subcontinental mantle material (Marroni et al., 2017).
Moreover, the record of high pressure rocks in the Sesia unit, reaching eclogite facies conditions (e.g. Regis et al., 2014), suggest subduction of part of the passive continental margin. More specifically, higher pressure conditions recorded in the former proximal part of the passive margin (Manzotti et al., 2014) indicates that underthrusting started in the middle of the margin, probably close to the necking zone (Fig. 6a) at ~90-95 Ma. The basal units of the Sesia-Dent Blanche nappes were also affected by intense ductile shearing associated with pre-Alpine deformation (Beltrando et al., 2010). Observations from the Alps are in agreement with our modelling scenario where the initiation of subduction is driven by a basal thrust at the base of the continental crust along the passive margin, subsequently crosscut and duplicated in the overall nappe system (Fig.   6). The original reactivation of the serpentinized mantle during contraction enabled the propagation of deformation toward the necking zone. Therefore, our study agrees that the necking domain localises subduction initiation, as previously inferred (Mohn et al., 2012), where the ductile lower crust is sufficiently thick to accommodate shortening (Fig. 6b). Our modelling also infers that a reactivation of the detachment at the base of the Sesia-Dent Blanche, Adriatic-derived extensional allochthon is not required to explain subduction initiation at the former Alpine passive margin (Manzotti et al., 2014).

Subduction initiation along the former Adriatic MRPM presently exposed in the Dinarides-Hellenides
Numerous studies have analysed the definition of tectonic units and their quantitative reconstruction in the Dinarides-Hellenides system of Central and SE Europe (e.g., Csontos & Vörös, 2004;Stampfli & Hochard, 2009). The most recent quantitative review  and reconstruction (van Hinsbergen et al., 2020) account for a separation between the European and Adriatic continent by the onset of the opening of the northern branch of the Neotethys (or Vardar) Ocean at ~250-245 Ma (Fig. 9b-c, e.g., Schmid et al., 2020). These times are when large volumes of Middle -Late Triassic intrusive and extrusive magmatism were documented along the Adriatic continental margin, from the Hellenides (e.g., Pe-Piper, 1998;Tsikouras et al., 2008) to the Dinarides (e.g., Pamic, 1984;Trubelja et al., 2004)

and the eastern Southern
Alps transitional area to the Alpine Tethys margin (e.g., De Min et al., 2020;Schmid et al., 2004a). The abundance of this basic to intermediate magmatism sourced from deep crustal or mantle levels (e.g., De Min et al., 2020;Pamic, 1984), associated with normal faulting and syn-kinematic deposition (e.g., Crisci et al., 1984;Robertson et al., 2009) infers partial melting during mantle upwelling and active rifting (e.g., Beltrán-Triviño et al., 2016). These observations combined with the Middle-Late Triassic age of calc-alkaline intermediate and mafic volcanism found in thrust ophiolitic mélanges (e.g. Pamic, 1984) testify the high magmatic budget along the Adriatic continental passive margin of the Dinarides-Hellenides (Fig. 10b). The subsequent closure of the northern Neotethys ocean by intra-oceanic subduction started during the late Early to Middle Jurassic, parallel to the ridge (Fig. 9a, ~160-180 Ma), inferred from dating metamorphic soles and mid-oceanic ridge (MOR) ophiolites, from occurrences of Late Triassic to Early Jurassic slices of scraped cherts and MOR basalts in thrust ophiolitic melanges and from the Middle-Late Jurassic age of the obducted ophiolitic sheet (Bortolotti et al., 2013;Mikes et al., 2008;Maffione et al., 2015;Ustaszewski et al., 2009). Similar to the Adriatic continental rifted margin facing the Alpine Tethys, spreading in the Neotethys ocean lasted for 50-70 Myr until subduction initiation (e.g. Schmid et al., 2004b, van Hinsbergen et al., 2020).
Plate reconstruction and paleo-magnetic studies suggest a convergence rate in the Neotethys Ocean of 1.5-2 cm yr -1 during 20 Myr (Capitanio & Goes, 2006;Handy et al., 2010;Le Breton et al., 2020van Hinsbergen et al., 2020. This ocean was ultimately closed by the onset of the latest Cretaceous continental collision and the formation of the Sava suture zone between Adriaticand European-derived continental units (e.g., Pamic, 1984;Schmid et al., 2020).
The comparison with our modelling study infers that the high magmatic budget of the Neotethys Adriatic continental rifted margin favoured the intraoceanic subduction (Fig. 7). This mode of subduction initiation is driven by the strong rheological coupling of the passive margin, inhibiting the continent-ward propagation of deformation (Fig. 7). The subduction initiation took place during 10-15 Myr (Bortolotti et al., 2013), which is in agreement with our modelling results (Figs. 3b). The age of oceanic relics found in ophiolitic melanges, combined with the age of obducted ophiolites, indicate that subduction initiated within the Early Jurassic part of the oceanic plate, although the occurrence of mid-ocean ridge ophiolites in contact with supra subduction zone ophiolites of the same age may also point to a subduction initiation close to the mid-oceanic ridge (Maffione et al., 2015). Although our models do not include a mid-oceanic ridge, they show that most shortening is initially localized in the place where the strength contrast is highest, i.e., at the ocean-continent transition, but migrates subsequently towards the oceanic domain, where the subduction zone forms at distances of less than 400 km from the margin (Fig. 5d).
Furthermore, our modelling predictions are comparable with previous observations in many other orogens containing ophiolitic belts, such as Anatolia (van Hinsbergen et al., 2016), Oman (Rioux et al., 2016;Saddiqi et al., 2006) or New-Caledonia (Lagabrielle et al., 2013). These observations have shown that the short time span (<10 Myr) between supra-subduction spreading centres and ophiolites emplacement points towards intra-oceanic subduction initiation close to the margin. These observations are also in agreement with our inferences showing that a coupled margin deforms significantly less at the ocean-continent transition and a subduction zone develops further into the ocean. We also acknowledge the potential role of plate re-organization and oblique convergence (Crameri et al., 2020;van Hinsbergen et al., 2020), or mantle rejuvenation by plume activity (Rolland et al., 2009), processes which are not accounted in our study.

Conclusions
Our thermo-mechanical models show that magma-poor passive margins with hyper-extended crust and exhumed serpentinized mantle lithosphere are favourable sites for induced subduction initiation. The serpentinized mantle facilitates strain localization and progressive development of a major shear zone that ultimately links up with zone of high strain that developed within the ductile lower continental crust in the necking zone. Shearing then propagates into the mantle lithosphere to initiate a subduction plate boundary. In contrast, passive margins intruded and underplated by magmatism favour intra-oceanic subduction initiation. These margins are rheologically coupled and stronger, which induces distributed deformation during the early shortening, while the subsequent deformation migrates within the oceanic lithosphere. These results indicate that the rheological and geometric configuration of the passive margin, inherited form the magmatic budget during continental rifting phase, is critical for strain localisation mechanisms during the early stages of deformation and shows that subduction initiation is mechanically feasible for Moho temperatures below 500 o C.
Our modelling results can be applied to the unique case of subduction initiation along two contrasting types of margins associated with the same Adriatic continental microplate. These results provide a mechanical explanation for subduction initiation at the ocean-continent transition in the western Alps and within the Neotethys ocean in the Dinarides-Hellenides orogenic systems, representing end-members of subduction initiation at magma-poor and magmarich margins, respectively.

Acknowledgments, Samples, and Data
The research project was funded by the European Union's EU Framework Programme for Research and Innovation Horizon 2020 "Subitop" under Grant Agreement No 674899. We are indebted to the editor and the reviewers for their valuable comments and suggestions, which have significantly improved the original manuscript.

Code availability
Researchers interested in using the numerical code should contact the authors. The code is available from the authors upon reasonable request.