# Buoyancy versus shear forces in building orogenic wedges

Lorenzo G Candioti, Thibault Duretz, Evangelos Moulas, Stefan M Schmalholz

► To cite this version:

Lorenzo G Candioti, Thibault Duretz, Evangelos Moulas, Stefan M Schmalholz. Buoyancy versus shear forces in building orogenic wedges. Solid Earth, 2021, 12 (8), pp.1749-1775. 10.5194/se-12-1749-2021. insu-03321729

HAL Id: insu-03321729

https://insu.hal.science/insu-03321729

Submitted on 18 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Buoyancy versus shear forces in building orogenic wedges

Lorenzo G. Candioti<fnref n="1" />, Thibault Duretz<fnref n="2" />, Evangelos Moulas<fnref n="3" />, and Stefan M. Schmalholz<fnref n="1" />

**Correspondence:** Lorenzo G. Candioti (lorenzo.candioti@unil.ch)

Received: 22 December 2020 – Discussion started: 23 February 2021

Revised: 9 June 2021 – Accepted: 17 June 2021 – Published: 10 August 2021

**Abstract.** The dynamics of growing collisional orogens are mainly controlled by buoyancy and shear forces. However, the relative importance of these forces, their temporal evolution and their impact on the tectonic style of orogenic wedges remain elusive. Here, we quantify buoyancy and shear forces during collisional orogeny and investigate their impact on orogenic wedge formation and exhumation of crustal rocks. We leverage two-dimensional petrological-thermomechanical numerical simulations of a long-term (ca. 170 Myr) lithosphere deformation cycle involving subsequent hyperextension, cooling, convergence, subduction and collision. Hyperextension generates a basin with exhumed continental mantle bounded by asymmetric passive margins. Before convergence, we replace the top few kilometres of the exhumed mantle with serpentinite to investigate its role during subduction and collision.

We study the impact of three parameters: (1) shear resistance, or strength, of serpentinites, controlling the strength of the evolving subduction interface; (2) strength of the continental upper crust; and (3) density structure of the subducted material. Densities are determined by linearized equations of state or by petrological-phase equilibria calculations. The three parameters control the evolution of the ratio of upward-directed buoyancy force to horizontal driving force, $F_B/F_D = Ar_F$, which controls the mode of orogenic wedge formation: $Ar_F \approx 0.5$ causes thrust-sheet-dominated wedges, $Ar_F \approx 0.75$ causes minor wedge formation due to relamination of subducted crust below the upper plate, and $Ar_F \approx 1$ causes buoyancy-flow- or diapir-dominated wedges involving exhumation of crustal material from great depth (> 80 km). Furthermore, employing phase equilibria density

models reduces the average topography of wedges by several kilometres.

We suggest that during the formation of the Pyrenees $Ar_F \lesssim 0.5$ due to the absence of high-grade metamorphic rocks, whereas for the Alps $Ar_F \approx 1$ during exhumation of high-grade rocks and $Ar_F \lesssim 0.5$ during the post-collisional stage. In the models, $F_D$ increases during wedge growth and subduction and eventually reaches magnitudes ($\approx 18 \text{ TN m}^{-1}$) which are required to initiate subduction. Such an increase in the horizontal force, required to continue driving subduction, might have “choked” the subduction of the European plate below the Adriatic one between 35 and 25 Ma and could have caused the reorganization of plate motion and subduction initiation of the Adriatic plate.

## 1 Introduction

The formation of collisional orogenic belts is an impressive manifestation of plate tectonics, and many studies have investigated the mechanisms causing mountain building in collisional settings (e.g. Malinverno and Ryan, 1986; Platt, 1986; Beaumont et al., 1996; Malavieille, 2010; Jaquet et al., 2018; Dal Zilio et al., 2020b). A popular geodynamic model explaining the formation of collisional orogens, such as the Western Alps, the Pyrenees or the Himalayas, is the wedge model (e.g. Chapple, 1978; Dahlen et al., 1984; Platt, 1986; Willett et al., 1993; Vanderhaeghe et al., 2003; Malavieille, 2010; Dal Zilio et al., 2020b). The first mechanical models of so-called critical wedges considered a frictional deformation (stresses are controlled by a specific yield criterion) and were

<footnote n="1">¹Institut des sciences de la Terre, Bâtiment Géopolis, Quartier UNIL-Mouline, Université de Lausanne, 1015 Lausanne (VD), Switzerland</footnote>

<footnote n="2">²Univ Rennes, CNRS, Géosciences Rennes, UMR 6118, 35000 Rennes, France</footnote>

<footnote n="3">³Institute of Geosciences & Mainz Institute of Multiscale Modeling (M³ODEL), Johannes-Gutenberg University, 55128 Mainz, Germany</footnote>
originally applied to accretionary wedges. The formation of such wedges has been extensively studied with both analogue and numerical models (e.g. Gutscher et al., 1998; Simpson, 2009; Malavieille, 2010; Graveleau et al., 2012; Ruh et al., 2012; Borderie et al., 2018; Dal Zilio et al., 2020b). Wedge models typically consider crustal deformation only and are driven by a kinematic boundary condition at the base of the crust, involving a rigid indenter, or backstop, which creates a kinematic singularity point at the base of the wedge. Such wedge models have also been used to study the formation of viscous fold nappes during fold and thrust belt evolution (e.g. Bauville and Schmalholz, 2015; Spitz et al., 2020) or the impact of surface processes on wedge formation (e.g. Willett, 1999).

Crustal wedge models have also been applied to entire collisional orogens and are frequently referred to as orogenic wedge models (e.g. Platt, 1986; Willett et al., 1993; Dal Zilio et al., 2020b). However, the geodynamic evolution of collisional orogens, such as the Alps or the Himalayas, typically involves the closure of oceanic domains and the subduction of oceanic and continental rocks before actual collision. Also, large regions of some of these orogens are characterized by exhumed high-pressure (> 1 GPa), and sometimes ultrahigh-pressure (> 2.7 GPa), rocks with peak temperatures ranging typically from 500 to 700°C (e.g. Lardeaux, 2014). Furthermore, tomographic images from the Western Alps and the Pyrenees (Zhao et al., 2015; Schmid et al., 2017; Teixell et al., 2018) indicate that the overriding lithospheric mantle is involved in the formation of orogenic wedges (see Fig. 1), implying that these orogenic wedges are rather lithospheric wedges (Nicolas et al., 1990) and not just crustal wedges. The above-mentioned pre-collisional subduction, the associated formation and exhumation of (ultra)high-pressure rocks, and the density contrast between the subducted crustal material and the surrounding mantle might significantly impact on the orogen dynamics (e.g. Beaumont et al., 1996; Burov et al., 2014; Butler et al., 2014; Sizova et al., 2014). Particularly, the subduction-related burial of crustal material to depths deeper than the isostatically balanced depth will cause (upward-directed) buoyancy forces, which act against the forces driving subduction and may assist rock exhumation.

From a mechanical point of view, the relatively slow tectonic deformation (no inertia) of the lithosphere is controlled by the balance of gravitational forces (acting everywhere inside a representative rock volume) and shear forces (or surface forces, acting on the surface of a representative rock volume) (e.g. Forsyth and Uyeda, 1975; Ramberg, 1981; Weijermars and Schmeling, 1986; Turcotte and Schubert, 2014; Gerya, 2019). In the absence of volumetric deformation, shear forces include forces acting normal or tangential to surfaces of the representative volume and induce shear deformation, such as pure or simple shear. Crustal wedge models typically consider shear forces and gravitational forces. The gravitational forces impact on frictional deformation via

a depth-increasing confining pressure and the gravitational adjustment of topographic gradients. Due to the kinematic singularity point that exists in most of the crustal wedge models, such models cannot predict deep subduction of continental crustal rocks. Consequently, (upward-directed) buoyancy forces resulting from density differences between subducted crust and surrounding mantle are not considered in these models. Here, we study the formation of orogenic wedges in a large-scale lithosphere (crust and mantle lithosphere) and upper mantle (asthenosphere and transition zone) framework, including subduction, to investigate the impact of shear and buoyancy forces on orogenic wedge formation.

We further aim to investigate orogenic wedge formation relevant to the Alpine orogeny. The formation of many collisional orogens, such as the Alps, is embedded in larger geodynamic cycles (Wilson, 1965; Wilson et al., 2019) including pre-orogenic rifting events. The Alpine orogen (see Fig. 1b) resulted from subduction of the Piemont–Liguria oceanic domain and collision of the European and Adriatic passive continental margins (e.g. Handy et al., 2010). Ophiolitic units found in the Alps mainly consist of serpentinized mantle material, indicating a Piemont–Liguria domain with exhumed and serpentinized mantle lithosphere (e.g. Manatschal and Müntener, 2009; McCarthy et al., 2020). For such Alpine-type orogens, structures inherited from continental margin formation presumably had a strong control on the dynamics of orogen formation (Chenin et al., 2017, 2019). Furthermore, progressive subduction of serpentinites and sediments may weaken the plate interface (Shreve and Cloos, 1986; Lamb and Davis, 2003; Behr and Becker, 2018) and impact on the balance between shear and buoyancy forces. Therefore, we consider models which involve a pre-collisional rifting phase and subduction involving serpentinites and sediments.

Here, we quantify the relative dominance and magnitude of buoyancy and shear forces during the formation of Alpine-type collisional orogens. We leverage two-dimensional (2D) petrological–thermomechanical numerical models of a long-term (ca. 170 Myr) extension–cooling–convergence cycle including subduction and continent–continent collision coupled to petrological-phase equilibria modelling (Gerya et al., 2004; Kaus et al., 2005; Yamato et al., 2007; Hetényi et al., 2007; Rummel et al., 2020). To quantify the relative importance of buoyancy and shear forces at convergent plate boundaries, we compare (1) simple or linearized equation of state (LEOS) density models to complex density (CD) models that include metamorphic reactions. The density structure of subducted continental crust and the surrounding mantle controls the magnitude of buoyancy forces. Increased (upward-directed) buoyancy forces eventually lead to slab detachment, vanishing of slab-pull forces (Duretz et al., 2012) and a rearrangement of forces throughout orogeny (Dal Zilio et al., 2020a). (2) We vary the strength of the upper crust and of serpentinites forming the subduction interface. The shear resistance of serpentinite lubricating the subduc-
<figure><img src="image_3.png" /></figure>

**Figure 1.** Simplified tectonic cross sections: (a) across the Central Pyrenees, modified after Teixell et al. (2018) and (b) across the Western Alps, modified after Zhao et al. (2015) and Schmid et al. (2017). Light violet represents the upper crust, dark violet represents the lower crust and green represents the mantle lithosphere. The region of mixed violet and green in (a) represents a unit of lower crust or mantle thrust slices. Black lines indicate unit boundaries and the grey line in (b) indicates the domain captured by most crustal wedge models (compare to Malavieille, 2010; Dal Zilio et al., 2020b, for example). Abbreviations: Penninic Front, PF; Briançonnais, BR; Dora-Maira, DM; and high-pressure regions, HP.

tion interface (Hirth and Guillot, 2013; Guillot et al., 2015, and references therein) may facilitate subduction. Also, coupling of the continental crust to the subducting mantle lithosphere impacts on the subduction and, therefore, orogen dynamics (Duretz and Gerya, 2013). We show that the relative dominance of shear or (upward-directed) buoyancy forces generates different modes, or styles, of orogenic wedges, some dominated by buoyancy-driven return flow and some by stacking of thrust sheets. We further analyse the importance of the shear resistance of serpentinites and the upper continental crust as well as rock density to make another step towards understanding the dynamics of Alpine-type orogens.

## 2 Numerical model

### 2.1 Mathematical model and numerical algorithm

As it is commonly done in continuum mechanics, the applied numerical algorithm solves the continuity and momentum equations coupled to conservation of energy expressed with respect to temperature. We consider incompressible visco-elasto-plastic materials that slowly flow (in-

ertia forces are negligible) under gravity and the influence of boundary tractions. Here, the term plastic refers to all rate-independent (instantaneous) irreversible deformation controlled by a Drucker-Prager yield function (Eq. A13). The applied pressure-dependent yield function is motivated by laboratory experiments (Byerlee, 1978) and is representative of rock deformation mechanisms at low temperatures (shear banding, faulting). The term viscous comprises all rate-dependent inelastic deformation, which is generally thermally activated. Material properties are advected using a marker-in-cell approach (Gerya and Yuen, 2003) involving initially 16 Lagrangian markers per Eulerian finite-difference grid cell. Up to 56 million numerical markers are advected in total. We do not apply any frictional or viscous strain softening. However, thermal softening, caused by shear heating, is active, resulting from the conservation of energy and temperature-dependent viscous flow stress. The algorithm has already been used to model deformation processes at various scales (Yamato et al., 2015; Duretz et al., 2016b; Yamato et al., 2019; Petri et al., 2019; Bessat et al., 2020) including upper mantle convection coupled to lithospheric-scale defor-
mation (Candioti et al., 2020). Appendix A provides a detailed description of the algorithm.

## 2.2 Model configuration

We here follow the modelling approach of Candioti et al. (2020) and we employ a 1600 km wide and 680 km deep (see Fig. 2) model domain. The global model resolution is 1 km × 1 km in horizontal *x* and in vertical *z* direction, respectively. A stress-free surface (Duretz et al., 2016a) is set initially at *z* = 0 km and, to allow for dynamic build-up of topography, the topmost +20 km of the model domain is left free of material. We impose constant material inflow and outflow velocities at the lateral boundaries, and the mechanical boundary at the bottom of the domain is free to slip. Further, we assume that no heat flows laterally out of the model domain, and the top and bottom temperature is kept constant at 15 and 1613 °C, respectively. Viscous deformation is modelled leveraging a combination of rheological flow laws, namely dislocation, diffusion and Peierls creep. Diffusion creep flow laws (Eq. A9) require knowledge of the grain size. Including an evolving grain size in lithospheric-scale models is usually achieved by employing additional parameterizations (e.g. paleowattmeters; Austin and Evans, 2007). As these parameterizations are commonly determined in laboratory experiments conducted at relatively fast deformation rates (compared to tectonic deformation), these parameterizations have to be extrapolated to natural conditions. This extrapolation introduces additional uncertainties. To keep the model simple, we here neglect grain size evolution and choose a constant average mantle grain size of 1 mm in Eq. (A9). Previous work has demonstrated that this choice, in combination with suitable flow law parameters, results in a viscosity structure that is consistent with geophysical constraints with respect to the convection dynamics of the mantle and the thermal thickness of the overlying, horizontal lithosphere (Candioti et al., 2020). Similar average grain sizes of the order of 1 mm in the upper mantle have been predicted by whole mantle convection models that include state-of-the-art grain size evolution models (Dannberg et al., 2017; Schierjott et al., 2020).

The effective density of the materials is either directly calculated via a linearized equation of state or pre-computed based on equilibrium phase-diagram sections for specific bulk-rock composition (see Appendix A). We consider surface processes such as erosion (constant erosion rate of 0.5 mm yr⁻¹ above *z* ≥ 2 km of topographic elevation) and sedimentation (instantaneous basin fill below *z* ≤ −5 km of topographic elevation, alternating calcite and mica rheology every 2 Myr). For more detailed justifications of model assumptions and explanations for the implementation of boundary conditions, the reader is referred to the study of Candioti et al. (2020). Model units include an initially 25 km thick upper crust with two vertical levels of 11 elliptical inclusions each, whose rheology is different from the upper

crustal matrix (weakened wet anorthite or Westerly granite rheology; see Table A1 for all material parameters). These elliptical inclusions represent structural and petrological units inherited from previous orogenic cycles and are motivated by field observations (e.g. Petri et al., 2019). Weak (wet quartzite) and strong (Maryland diabase) inclusions may represent metasedimentary units and mafic intrusions, respectively. Below the upper crust an 8 km thick lower crust (wet anorthite) is employed. The initial Moho is set to *z* = −33 km and the mantle lithosphere (dry olivine) extends down to *z* = −120 km. In total, 12 weak (wet olivine) elliptical inclusions over two vertical levels are included in the mantle lithosphere. These heterogeneities are motivated by geophysical observations (e.g. Barnhoorn et al., 2010) and represent weaknesses owing to spatial lithological variations. Elongated or elliptical inclusions allow for the generation of asymmetric margin geometries (Duretz et al., 2016b; Petri et al., 2019). We include the asthenosphere and transition zone (grouped and termed upper mantle hereafter, dry olivine rheology) down to *z* = −660 km. The distribution and emplacement of the inclusions is described in detail in Appendix B. The difference between mantle lithosphere and upper mantle is due to temperature and pressure only; i.e. all material parameters are the same.

## 2.3 Simulations

Four types of simulations are performed: (1) the reference model (REF) generates a long-term (> 160 Myr) geodynamic cycle of subsequent extension-cooling-convergence in a single and continuous simulation. During a 50 Myr rift-ing period, an ≈ 360 km wide basin is formed. This basin is floored by exhumed continental mantle and bounded by two hyper-extended magma-poor rifted margins. We apply an absolute extension velocity of 1 cm yr⁻¹ and we assume an ultra-slow to slow spreading rift system with minor melt production and therefore neglect decompressional melting of the peridotites. No far-field plate velocity is applied to the system for the following 60 Myr. This choice of applied deformation velocity and duration of the periods leads to margin geometries comparable to reconstructions from the ancient Alpine Tethys margin system and Piemont–Liguria ocean (Le Breton et al., 2021). In all models presented here, the subsequent convergence is driven by a kinematic boundary velocity of 1.5 cm yr⁻¹ (absolute value). The convergence phase involves subduction initiation, closure of the basin floored by exhumed mantle, continental collision and orogenic wedge formation. Model REF simulates a scenario without serpentinization of the exhumed mantle peridotite. Except for the mantle lithology that utilizes a complex density (CD) model, the density of all other model lithologies is calculated using a LEOS (Eq. A3).

The geodynamic evolution of REF at the end of the cooling stage (109 Myr) serves as a self-consistently generated initial configuration for three additional types of simulations.
<figure><img src="image_4.png" /></figure>

**Figure 2.** Model configuration and boundary conditions. (a) Profile of horizontal velocity for material inflow and outflow along the western model boundary. The blue line indicates the profile for the extension, the purple line indicates the profile for the cooling and the yellow line indicates the profile for the convergence. (b) Entire model domain. Crustal phases: yellow represents the crustal matrix, dark orange represents strong elliptical inclusions and light orange represents weak elliptical inclusions. White to red shows the effective viscosity field as calculated by the numerical algorithm. All colour maps used to visualize physical fields in this study are provided by Crameri (2018). Black line shows the vertical temperature profile. (c) Enlargement of the central region of the domain. (d) Same profile as shown in (a), but along the eastern model boundary.

(2) We parameterize a serpentinization front propagating through the topmost layer of exhumed peridotites to test the impact of serpentinite strength on the convergent deformation. In this parametrization, we replace the dry olivine flow law parameters of the exhumed peridotites in the basin above $z = -9$ km with parameters for antigorite rheology (Hilairet et al., 2007). This approach results in an effective average thickness of the serpentinite layer of ca. 5.5 km. In order to investigate the effects of variable serpentinite strength, we gradually increase the prefactor in front of the dislocation creep flow law ($\zeta$ in Eq. A8) for antigorite (Hilairet et al., 2007). We here report the modelling results for prefactors 1 and 18, because these models are end-member types for subduction and orogen dynamics. The background upper and lower crust are feldspar-dominated (see Table A1 for details). As in REF, density is computed with a LEOS except for the mantle where a CD model is employed. (3) We consider serpentinization and investigate the impact of (upward-directed) buoyancy forces on the collisional stage employing a CD model for all lithologies. Type (4) is identical to type (3), except that compared to a feldspar-dominated (wet anorthite) rheology a relatively weaker quartz-dominated rheology (Westerly granite) is employed for the upper crust. A summary of all simulations is given in Table 1. The applied material parameters are given in Table A1, and the bulk rock compositions as well as the solution models used for phase diagram calculations are given in Table 2.

## 3 Results

We first describe the results of the reference model, REF, for the entire geodynamic extension-cooling-convergence cycle. We then describe the impact of serpentinite shear resistance on the evolution of the models during convergence. Finally, we report results of linearized and complex density models as well as the impact of upper crustal shear resistance. We focus on the continent-continent collision stage of the individual models. Stages of basin closure, serpentinite channel formation and the subduction dynamics are documented in the video supplement of models REF (Candioti, 2020a) and GC1 (Candioti, 2020b) but are not described in further detail here.

### 3.1 Reference model – REF

Figure 4 shows the geodynamic evolution of the reference model (REF). Crustal break-up occurs at $\approx 19$ Myr in model history (Fig. 4a). At this stage, hot mantle material rises in the horizontal centre and diverges laterally below the plates (see arrows in Fig. 4a). The 1300 $^{\circ}$C isotherm remains at a depth of $\approx 20–30$ km below the rift centre (Fig. 4d). Two asymmetric continental margins have formed. The necking zone of the left margin is $\approx 20$ km wide. In total, the width of the left margin is ca. 80–90 km. The right margin is in total ca. 140–160 km wide including an $\approx 60–80$ km wide necking zone. At the end of the cooling period (109 Myr, Fig. 4b), an $\approx 360$ km wide basin has opened, which is floored by exhumed mantle material. Convection in the upper mantle stabilizes the thermal and mechanical thickness of the litho-
<table><caption>**Table 1.** Summary of all simulations presented in this study. Bold entries highlight the differences compared to the reference model REF.</caption><thead><tr><th>Parameter</th><th>Unit</th><th>REF</th><th>AS1</th><th>AS18</th><th>AC1</th><th>AC18</th><th>GC1</th><th>GC18</th></tr></thead><tbody><tr><td>Crustal matrix*</td><td>–</td><td>Anorthite</td><td>Anorthite</td><td>Anorthite</td><td>Anorthite</td><td>Anorthite</td><td><strong>Granite</strong></td><td><strong>Granite</strong></td></tr><tr><td>Complex density (CD) model</td><td>–</td><td>No</td><td>No</td><td>No</td><td><strong>Yes</strong></td><td><strong>Yes</strong></td><td><strong>Yes</strong></td><td><strong>Yes</strong></td></tr><tr><td>Serpentinization front</td><td>–</td><td>No</td><td><strong>Yes</strong></td><td><strong>Yes</strong></td><td><strong>Yes</strong></td><td><strong>Yes</strong></td><td><strong>Yes</strong></td><td><strong>Yes</strong></td></tr><tr><td>Serpentinite pre-factor (ζ)</td><td>–</td><td>–</td><td><strong>1</strong></td><td><strong>18</strong></td><td><strong>1</strong></td><td><strong>18</strong></td><td><strong>1</strong></td><td><strong>18</strong></td></tr><tr><td>Starting time</td><td>Myr</td><td>0</td><td><strong>109</strong></td><td><strong>109</strong></td><td><strong>109</strong></td><td><strong>109</strong></td><td><strong>109</strong></td><td><strong>109</strong></td></tr></tbody></table>

sphere to ca. 90–100 km (no velocity glyphs in Fig. 4b and e; see also the discussion in Candioti et al. (2020) for more detail). At 112 Myr subduction is initiated at the transition from the proximal to the distal continental margin in the right model side (see Sutra et al., 2013, for nomenclature). Partly, the necking zone of this margin has been subducted to ca. 40 km depth (Fig. 4c and f). At 144 Myr in model history, the continental crust of the left margin (blue colours in Fig. 5a and d) is subducted to ca. 60 km depth. Sediments originally deposited in the trench have been incorporated in between the subducting and the overriding plate. With ongoing convergence, the upper crust of the subducting plate is buried to a maximum depth of ca. 120 km. At 155 Myr, upper crustal units are sheared off the subducting plate at the transition to the lower crust and begin to form thrust sheets (see deflected isotherms at $x \approx -80$ km and $z \approx -30$ km in Fig. 5e). Major parts of continental crust are sheared off the subducting plate before entering the subduction zone. An orogenic wedge with several thrust sheets has formed at 164 Myr (Fig. 5f). At this stage, the lower part of the subducting slab detaches and the deeper subducted upper crustal material flows upward. This return flow is limited to depths between ca. 60 and 80 km (see arrows in Fig. 5f).

## 3.2 Impact of low (AS1) and high (AS18) serpentinite shear resistance

As in REF, subduction is initiated at the transition from the proximal to the distal right continental margin in AS1 and AS18. A large volume of serpentinite material is sheared off the subducting plate in AS1 and AS18 (Fig. 6a and e). In AS1, the serpentinites form a coherent, inclined channel down to ca. 100–120 km depth (Fig. 6a). At 148 Myr, the upper continental crust of the left margin (blue colours in Fig. 6) has been subducted to ca. 90–100 km depth (Fig. 6a). Between 154 and 162 Myr, the subducted upper continental crust returns to the surface through the weak serpentinite channel in AS1. An $\approx 50$ km wide block of the overriding plate’s upper crust is separated from the right margin by subducted crustal material that is being exhumed to the surface. At 164 Myr, large parts of the subducted continental crust have been exhumed to the surface (Fig. 6d). Units containing weak inclusions are internally deformed during the ex-

humation. Strong units are folded but exhumed coherently (Fig. 6d). The serpentinite material is surrounding the exhumed upper crustal material and the separated continental upper crustal block of the overriding plate. In AS18, the serpentinite channel is intersected by the two colliding plates at $z \approx -30$ km (Fig. 6e). In contrast to AS1, the continental upper crust is partly sheared off the subducting plate at $z \approx -25$ km in AS18 (Fig. 6e). Maximum burial depth of upper continental crust in this model is ca. 80–90 km. Exhumation of subducted crustal material is limited to $z \approx -70$ km (Fig. 6f and g). Instead, the upper continental crust is wedged at shallower depths, similar to REF (see deflected isotherms at $x \approx -125$ km in Fig. 6g). Thrust sheets propagate along the upper continental crust of the subducting plate (Fig. 6h). The serpentinized mantle material forms a coherent unit emplaced in between the wedged continental upper crust of the subducting plate and the relatively undeformed continental upper crust of the overriding plate. Dimensions of the evolved orogenic wedge in AS18 are similar to the wedge dimensions evolved in REF (compare Fig. 5f to Fig. 6h).

## 3.3 Impact of complex density models – models AC1 and AC18

Density models that include predicted metamorphic assemblages (CD models) lead to more variable density fields and larger density changes as a function of pressure and temperature compared to LEOS models (compare Fig. 3g to f). Similar to AS1, a serpentinite channel (low viscosity region in Fig. 7a) forms after subduction initiation in AC1 (Fig. 8a). The continental upper crust of the subducting plate (blue colours) is subducted to ca. 150 km depth, which is deeper compared to AS1 (compare Fig. 6b to Fig. 8b and c). In contrast to AS1, exhumation of the subducted upper continental crust has not occurred until 163 Myr (compare Fig. 6c to Fig. 8c) in model history. At 170 Myr, the subducted continental crust in AC1 has been exhumed to the surface along the weak subduction interface (Fig. 8d). In contrast to AS1, the continental upper crust of the overriding plate has not been separated by the returning continental crust of the subducting plate in AC1 (compare Fig. 6d to Fig. 8d). Figure 9a and b show the density field of models AS1 and AC1, respectively, computed by the algorithm at 170 Myr in model

<footnote>* All material parameters are listed in Table A1.</footnote>
<table><caption>**Table 2.** Bulk rock composition for phase equilibrium calculations.</caption><thead><tr><th>Oxides [wt %]</th><th>Pelite (avg.)<sup>1,4</sup></th><th>Rhyolite<sup>1,4</sup></th><th>Andesite<sup>1,4</sup></th><th>MORB<sup>1,4</sup></th><th>Hydr. peridotite<sup>1,4</sup></th><th>Serpentinite<sup>2,4</sup></th><th>Bulk DMM<sup>3,5</sup></th></tr></thead><tbody><tr><td>SiO<sub>2</sub></td><td>61.500</td><td>72.800</td><td>57.900</td><td>49.200</td><td>44.710</td><td>44.210</td><td>44.710</td></tr><tr><td>Al<sub>2</sub>O<sub>3</sub></td><td>18.600</td><td>13.300</td><td>17.000</td><td>16.100</td><td>4.160</td><td>3.130</td><td>3.980</td></tr><tr><td>FeO</td><td>10.000</td><td>2.440</td><td>6.980</td><td>10.220</td><td>8.070</td><td>8.898</td><td>8.180</td></tr><tr><td>MgO</td><td>3.810</td><td>0.390</td><td>3.330</td><td>6.440</td><td>39.200</td><td>39.240</td><td>38.730</td></tr><tr><td>CaO</td><td>—</td><td>—</td><td>6.790</td><td>10.500</td><td>2.420</td><td>3.060</td><td>3.170</td></tr><tr><td>Na<sub>2</sub>O</td><td>1.460</td><td>3.550</td><td>3.480</td><td>3.010</td><td>0.220</td><td>—</td><td>0.130</td></tr><tr><td>K<sub>2</sub>O</td><td>3.020</td><td>4.300</td><td>1.620</td><td>1.100</td><td>—</td><td>—</td><td>—</td></tr><tr><td>H<sub>2</sub>O</td><td>sat</td><td>sat</td><td>sat</td><td>sat</td><td>sat</td><td>sat</td><td>sat</td></tr><tr><td colspan="8">Solution models</td></tr><tr><td>Op<sub>x</sub>(HP)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Gt(GCT)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>feldspar</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Chl(HP)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Sp(HP)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>O(HP)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>St<sub>lp</sub>(M)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Carp</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Sud</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Bio(TCC)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>St(HP)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Ctd(HP)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Pheng(HP)</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td><td>—</td></tr><tr><td>hCrd</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Omph</td><td>—</td><td>—</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>GlTrTsPg</td><td>—</td><td>—</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Pu(M)</td><td>—</td><td>—</td><td>+</td><td>+</td><td>+</td><td>—</td><td>—</td></tr><tr><td>Act(M)</td><td>—</td><td>—</td><td>+</td><td>+</td><td>+</td><td>+</td><td>—</td></tr><tr><td>T</td><td>—</td><td>—</td><td>+</td><td>—</td><td>+</td><td>+</td><td>—</td></tr><tr><td>A-phase</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Chum</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td><td>+</td><td>—</td></tr><tr><td>B</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Wus</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Fperh</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td><td>+</td><td>—</td></tr><tr><td>Atg(PN)</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td><td>+</td><td>—</td></tr><tr><td>C2/c</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Wus</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Pv</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Pl</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Sp</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>O</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Wad</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Ring</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Op<sub>x</sub></td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Cpx</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Aki</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Gt_maj</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>Ppv</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr><tr><td>CF</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>+</td></tr></tbody></table>

Modified bulk rock after ¹ Winter (2013), ² Pelletier et al. (2008), ³ Workman and Hart (2005). We assume water saturation in all calculations. Crosses denote solution models used for given lithologies. Thermodynamic databases used: ⁴ Holland and Powell (1998) updated in 2002 and ⁵ Stixrude and Lithgow-Bertelloni (2011) for depleted MORB mantle (DMM). Note that only the relative abundance of elements is important for our phase diagram calculations. Details on the solution models can be found in the solution_model.dat data file in Perple_X.
history. Compared to AS1, in AC1 the density of the ma-
terial in the wedge is much more variable. The density of
the subducted crustal units in AC1 is up to ca. 100 kg m⁻³
higher than in AS1. The resulting buoyancy contrasts are
much smaller, leading to significantly less uplift of topogra-
phy in AC1 compared to AS1. In AS1, the maximum eleva-
tion of topography exceeds 10 km; in AC1 the maximum el-
evation is locally ≈ 7 km, but high topographies are on aver-
age ≈ 5 km. The width of the wedge is ≈ 350 and ≈ 275 km
in AS1 and AC1, respectively. Compared to AS18, only a mi-
nor volume of serpentinite material has been sheared off the
subducting slab in AC18 (compare Fig. 6e to Fig. 8e). The
majority of serpentinite material is coupled to the slab and
subducted into the upper mantle (Figs. 8f and 7c). Subduc-
tion of continental upper crust reaches a depth of ca. 180 km
after 163 Myr (Fig. 8g). In the absence of a weak, serpen-
tinized subduction interface, the resistance of the subduction
channel is increased in AC18 (compare Fig. 7a to c). The
deeply subducted continental upper crust cannot be exhumed
along the strong subduction interface and breaks through the
less resistant mantle wedge (Fig. 8h), relaminating below the
overriding plate (similar to the models presented in Currie
et al., 2007, and Li and Gerya, 2009).

3.4 Impact of crustal strength – models GC1 and GC18

Similar to AS1 and AC1, in GC1 the serpentinites form a
channel along the subduction interface down to ca. 120 km
depth (Figs. 10a and 7b). Continental upper crust of the sub-
ducting plate is buried to z ≈ −120 km (Fig. 10a and b). At
c.a. 153 Myr, the subducting slab detaches at a depth of ca.
400 km (see Candioti, 2020b). Between 156 and 162 Myr,
the subducted continental upper crust flows back to the sur-
face along the weak subduction interface breaking through
the continental upper crust of the overriding plate (Fig. 10c).
The serpentinite material is surrounding the exhumed crustal
material (Fig. 10d) within the orogenic wedge. Exhumation
of upper continental crust from z ≈ −80 km to z ≈ −30 km
occurs between 162 and 164 Myr. At ca. 165 Myr, a sec-
ond slab detachment occurs at z ≈ −250 km (see Candioti,
2020b). In contrast to AC18, less crustal volume is subducted
in GC18 (Fig. 10e). Instead, the continental upper crust is
largely sheared off the subducting lithosphere at z ≈ −40 km
(Fig. 10f), and several thrust sheets form (Fig. 10f). Com-
pared to GC1, the relatively stronger serpentinite in GC18
is largely subducted (Fig. 7d) and sediments originally de-
posited in the trench are incorporated into a growing orogenic
thrust wedge (Fig. 10f–h).

**3.5 Evolution of buoyancy and shear forces**

In order to investigate the relative impact of buoyancy and
shear forces on collision, we quantify the temporal evolu-
tion of the horizontal driving force per unit length (Fig. 11b,
$F_D$ hereafter; see Appendix C) and of the (upward-directed)

buoyancy forces of subducted crustal material (Fig. 11c, $F_B$ hereafter; see Appendix C). Buoyancy forces calculated here represent an upward-directed (positive) buoyancy force of the subducted crust acting against further subduction. The ratio $F_B/F_D$ (Fig. 11a, $Ar_F$, hereafter) is a measure for the relative dominance of buoyancy or shear forces driving the deformation and exhumation within the orogenic wedge (see Appendix C). Points 1–5 in Fig. 11a–e represent important stages of model GC1 (thick dashed turquoise blue line) and are representative for the evolution of all presented models.

Subduction initiation occurs for $F_D \approx 22 \text{ TNm}^{-1}$ in models of feldspar dominated continental upper crust (AS1, AS18, AC1 and AC18) and for $F_D \approx 18 \text{ TNm}^{-1}$ in models of quartz dominated upper crust (GC1, GC18, 1 and grey area in Fig. 11b). Magnitudes of $F_D$ decrease after subduction initiation to ca. 10–15 TNm$^{-1}$ and remain, on average, relatively constant during basin closure (between 1 and 2 in Fig. 11). Until basin closure, $Ar_F \approx 0.25$ and remains relatively constant in all models. Magnitudes of $F_B$ increase with an increasing volume of crustal material involved in subduction (see Fig. 11c and d) between 140 and 153 Myr. In response, magnitudes of $Ar_F \approx 1$ are reached in AS1, AC1 and GC1. Between ca. 153 and 156 Myr, shallower buried crustal material is exhumed in GC1 (3 in Fig. 11; see also Fig. 10a) coinciding with a deep slab detachment (see Candioti, 2020b). At ca. 162 Myr, the deeper subducted units are exhumed (see also Fig. 10c). This exhumation event is followed by a shallow slab detachment leading to a rapid increase in magnitude of $F_D$ at ca. 165 Myr. The buoyancy pull of the subducting slab is lost, and a larger horizontal driving force is needed to overcome the (upward-directed) buoyancy push of subducted crust and to continue subduction with the prescribed kinematic boundary velocity. While the buoyancy forces remain relatively constant, the increase in $F_D$ decreases the magnitude of $Ar_F$ (between 4 and 5 in Fig. 11a–c). The maximum burial depth of continental crust is reached in AC1 and AC18 (feldspar-dominated upper crust) and is ca. 200 km. In REF, GC1, GC18, AS1 and AS18 the maximum burial depth varies between 120–150 km. During the formation of the collisional orogen, the magnitude of $F_D$ increases and eventually exceeds the magnitude that was necessary for subduction initiation (grey area in Fig. 11b). Magnitudes of $F_D$ and $F_B$ are on average higher in models of feldspar-dominated continental upper crust compared to models of quartz-dominated upper crust. The increased shear resistance in these feldspar-dominated models allows for deeper subduction of continental upper crust when the serpentinite is weak. High values for $F_D$ in AS18 and REF are explained by a more resistant subduction interface caused by either an increased shear resistance of the serpentinite material, or absence of serpentinization.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Density structures employed in the models presented here. (a–e, g–h) Phase diagram density fields calculated with Perple_X for the bulk rock compositions given in Table 2. (f, i) Density fields calculated using a linearized equation of state (Eq. A3) for an upper crustal reference density of 2800 kg m⁻³ and a serpentinite reference density of 2585 kg m⁻³, respectively.</figcaption></figure>

## 3.6 Different modes of orogenic wedge formation

The geometry, kinematics and dynamics of the evolving collision zone and orogenic wedge varies significantly in the presented models (Fig. 12). We refer here to these variations as different modes of orogenic wedge formation. These varying modes depend mainly on the shear resistance of upper crust and serpentinite, and its buoyancy contrast to the surrounding mantle material. The magnitude of ArF is used to classify the different modes of orogenic wedge formation.

For low shear resistance of serpentinites, generating a weak subduction interface, ArF ≈ 1, indicating equal importance of buoyancy and shear forces building the orogenic wedge. Two modes are observed: (1) a diapir-like mode of exhumation due to either (i) increased buoyancy contrast (LEOS models) and high shear resistance of the upper crust (AS1; see Fig. 12a) or (ii) a decreased buoyancy contrast (CD models) and low shear resistance of the upper crust (GC1; see Fig. 12c); (2) a channel flow mode of exhumation occurs for a decreased buoyancy contrast (CD models) and increased shear resistance of the upper crust (AC1; see Fig. 12b).

In models of high serpentinite shear resistance, generating a strong subduction interface, two different modes are also observed: (3) a thrust-diapir mode of deformation for either (i) an increased buoyancy contrast and high shear resistance of the upper crust or (ii) a decreased buoyancy contrast and low shear resistance of upper crust. In both cases ArF ≈ 0.5, indicating that shear forces dominate the building of the orogenic wedge (Fig. 12d and f).

(4) A relamination mode due to relamination of deeply subducted continental crust is observed for a decreased buoyancy contrast and high shear resistance of the upper crust (ArF ≈ 0.75; see Fig. 12e).

# 4 Discussion

## 4.1 Buoyancy vs. shear forces controlling modes of orogenic wedge formation

The shear resistance of (i) serpentinites and (ii) the upper continental crust directly impacts the shear forces (compare GC1 to AC1 in Fig. 11b). Determining the shear resistance, or effective viscosity, of the lithosphere deforming under geological timescales remains challenging (Burov and Watts, 2006). Rock deformation experiments are performed at deformation rates that are many orders of magnitude higher than tectonic deformation rates. Best-fitting curves have to be extrapolated to natural conditions, which introduces large uncertainties on the actual strength of rocks deforming under natural conditions (Mancktelow and Pennacchioni, 2010; Idrissi et al., 2020). Serpentinite plays a crucial role in subduction zones and, ultimately, in the formation of Alpine-type collisional orogens (Hess, 1955; Raleigh and Paterson, 1965; Hacker et al., 2003; McCarthy et al., 2020). Despite its importance, the rheology of serpentinite at lithospheric-scale pressure and temperature conditions remains elusive (David et al., 2018; Hirauchi et al., 2020, and references therein).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Evolution of the reference run REF during rifting, cooling and subduction initiation. Panels (a)–(c) show the entire domain and (d)–(f) show an enlargement. White to dark red indicates the effective viscosity of the mantle material calculated by the algorithm, yellow to orange indicates the crustal matrix and the weak and strong inclusions of the overriding plate, light to dark blue indicates the corresponding crustal units of the subducting plate. The material parameters used for the upper crustal phase of the overriding and the subducting plates are identical. Green indicates the lower crust of both plates and salmon and brown indicates the sedimentary units. White lines show different isotherms and the arrows indicate the velocity field calculated by the numerical algorithm.</figcaption></figure>

Several deformation mechanisms for serpentinite material, often based on experiments with antigorite, have been discussed in the literature, including dislocation creep (Hilairet et al., 2007), semi-brittle or plastic deformation behaviour (Chernak and Hirth, 2010; Hirth and Guillot, 2013), grain boundary sliding (Idrissi et al., 2020), and sliding on shear cracks (Hansen et al., 2020). Numerical models are useful to test different end-member rock rheologies and to investigate the impact of rock strength on subduction and formation of collisional orogens. To first order, weak serpentinite material may indeed form a subduction channel and lubricate the subduction interface. The subduction channel in model AC1 has formed self-consistently, because subduction was initiated without a prescribed major weak zone in the lithosphere. Deep subduction and exhumation to the surface of crustal material is feasible in this kind of model. In contrast, strong serpentinite material leads to detaching and thrusting of crustal material already escaping subduction at shallow depths. In natural settings, the effective strength of serpentinites may vary along the upper regions of the subducting plate parallel and/or orthogonal (along-trench) to the subduction direction, for example due to varying degrees of

serpentinization. Such spatial strength variation may cause temporal and/or along-trench alternations between channel mode and thrust-diapir mode when serpentinites of different strength enter the subduction zone. The shear resistance of the upper crust resulting from the rheological flow laws employed here, wet anorthite (Rybacki and Dresen, 2004) and Westerly granite (Hansen and Carter, 1983), is similar and neither extremely low nor extremely high (see also Fig. 1 in Bürgmann and Dresen, 2008). However, employing these two different flow laws for the upper crust changes the mode of orogenic wedge formation significantly. Similar to the strength of serpentinites, in natural settings the effective strength of the upper crust may also vary and cause temporal and along-trench variations of the orogenic wedge modes during the evolution of continental collision.

Magnitudes of buoyancy forces in our models are enhanced by large density contrasts between the subducted material and the surrounding mantle. We tested end-member models of simple and complex density calculations. The pre-computed density tables are based on calculated equilibrium phase diagram sections (Fig. 3). In our calculations we have assumed H₂O saturation and therefore our system is open
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Evolution of the reference run REF after closure of the marine basin and continent–continent collision. White to dark red indicates the effective viscosity of the mantle material calculated by the algorithm, yellow to orange indicates the crustal matrix and the weak and strong inclusions of the overriding plate, light to dark blue indicates the corresponding crustal units of the subducting plate. The material parameters used for the upper crustal phases of the overriding and the subducting plates are identical. Green indicates the lower crust of both plates and salmon and brown indicates the sedimentary units. White lines are different isotherms and the arrows indicate the velocity field.</figcaption></figure>

with respect to its H₂O content, but closed with respect to other elements. This is a valid approximation for the kilometre scale that we are considering. However, the preservation of high-grade metamorphic rocks at the near-surface environment indicates that the crustal metamorphic rocks are not always at thermodynamic equilibrium as is assumed here. The rate of metamorphic re-equilibration is strongly affected by temperature and by the availability of fluids (Rubie, 1986; Austrheim, 1987; Malvoisin et al., 2020). Coupling mineral-scale phase equilibria modelling to large-scale geodynamic models remains challenging. Although coupling of petrological and thermo-mechanical modelling via CD models as presented here is simplified, this approach has proven useful to explain observations that cannot be predicted by the commonly used LEOS models. These observations include (i) varying sediment thickness accumulated during basin subsidence (Kaus et al., 2005), (ii) evolution of subducting slab dynamics (van Hunen et al., 2001; Toussaint et al., 2004) and (iii) the exhumation of (U)HP rocks (Yamato et al., 2007; Warren et al., 2008). In addition, our results demonstrate that CD models avoid unrealistically high topographic elevation during orogen formation.

Varying mechanical strength of (i) serpentinites and (ii) the upper crust in combination with varying buoyancy contrasts changes the mode of orogenic wedge formation modelled here (Fig. 12). Hence, it is essential to (1) further refine our knowledge on the mechanical strength of crustal rocks and serpentinites under natural deformation conditions and (2) to account for realistic density structures including metamorphic reactions in numerical models of collisional orogen formation.

## 4.2 Subduction initiation

We follow here the modelling approach by Candioti et al. (2020). The margin geometry and thermal structure, prior to convergence, is generated during a modelled rifting and cooling period. This way, the generated passive margin, the marine basin system and the subsequent convergence are modelled in an internally consistent manner. Subduction initiation is horizontally forced (Stern, 2004; Stern and Gerya, 2018; Crameri et al., 2020), and a major lithospheric shear zone forms around the transition from the distal to the proximal margin (see also the discussion in Candioti et al., 2020). The ad hoc parameterized layer of serpentinite is not relevant
<figure><img src="image_8.png" /></figure>

**Figure 6.** Convergence and collisional stage of AS1 and AS18 (low and high shear resistance of serpentinite, LEOS density model, feldspar-dominated upper crust). White to dark red indicates the effective viscosity of the mantle material calculated by the algorithm, yellow to orange indicates the crustal matrix and the weak and strong inclusions of the overriding plate, light to dark blue indicates the corresponding crustal units of the subducting plate. The material parameters used for the upper crustal phases of the overriding and the subducting plates are identical. Green indicates the lower crust of both plates and salmon and brown indicates the sedimentary units. White lines are different isotherms and the arrows indicate the velocity field calculated by the numerical algorithm.

for subduction initiation, as subduction is initiated also in the reference model without serpentinite (see Fig. 4f). Instead, geometrical focusing of stresses below the margin together with thermal softening and a temperature-dependent viscosity leads to spontaneous formation of a shear zone transecting the lithosphere (Thielmann and Kaus, 2012; Jaquet and Schmalholz, 2018; Kiss et al., 2020; Candioti et al., 2020; Auzemery et al., 2020). In our models, magnitudes of $F_D$ between 18 and $22 \text{ TN m}^{-1}$ are necessary to initiate subduction for quartz- and feldspar-dominated upper crust, respectively. These magnitudes are significantly lower compared to $F_D \approx 37 \text{ TN m}^{-1}$ obtained by Kiss et al. (2020). They studied subduction initiation at an idealized, ad hoc constructed passive margin without mechanical heterogeneities in the form of a multi-layer or elliptical geometry (see also Candioti et al., 2020, for a detailed comparison). The Peierls flow law

parameters describing the rheology of olivine used here have been elaborated by Goetze and Evans (1979). Recent studies suggest that the olivine strength resulting from this parameterization is likely overestimated (Idrissi et al., 2016). If true, stresses in the mantle lithosphere would be lower than predicted by our models, which could further reduce the magnitude of shear forces necessary for subduction initiation. However, the minimum value of $F_D$ required for subduction initiation in our models should be determined in future studies.

## 4.3 Potential applications to natural collisional orogens

### 4.3.1 Implications for the Pyrenean orogeny

The opening of the Bay of Biscay (110–105 Ma, Henry et al., 1998) lead to the formation of a narrow marine basin which
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Plate-interface and upper crustal-matrix strength for models AC1, AC18, GC1 and GC18. White to dark red is the effective viscosity field calculated by the numerical algorithm. (a, b) Weak interface (AC1 and GC1), (c, d) strong interface (AC18, GC18). Grey dashed lines are several isotherms, cyan solid lines are material phase boundaries. Abbreviations: UC = upper crust, LC = lower crust and S = serpentinites. Low viscosity of serpentinites indicates weak plate interface (a, b). For high viscosity the serpentinites are coupled to the lower plate and subducted into the mantle (c, d). For a wet anorthite upper crustal matrix (a, c) the brittle-ductile transition (depth of plastic shear bands) in the upper crust is at greater depth compared to a Westerly granite upper crustal matrix (b, d).</figcaption></figure>

was floored by exhumed mantle in the present-day Pyrenean domain. Convergence between the Iberian and European plates initiated at ca. 85 Ma and culminated in the formation of the Pyrenean collisional orogen in the Eocene (Jammes et al., 2009). Orogenic wedge formation in the Pyrenees involved mainly the upper 20 km of the upper crust (see Fig. 1a; Teixell et al., 2018, and references therein). The absence of subduction-related high-pressure metamorphism (Muñoz, 1992) indicates that no significant volumes of upper crust have been subducted to great depth (> 50 km). Instead, the majority of the upper crust has been presumably sheared off the subducting lower crust and formed thrust sheets at mid to upper crustal level (Muñoz, 1992; Teixell et al., 2018, and references therein). Rifting-related inheritances are likely important for the formation of crustal thrust sheets during the Pyrenean orogeny (Jammes and Huismans, 2012; Jammes et al., 2014; Erdős et al., 2014). In our models, a ca. 360–400 km wide basin (see for example Fig. 4b) with exhumed serpentinized mantle is generated during a rifting period modelled prior to convergence. The model geometry before the onset of convergence is, thus, not directly applicable to the pre-orogen geodynamic setting in the Pyrenees, because the basin resulting from the opening of the Bay of Biscay was most likely considerably narrower (Jammes et al., 2009; McCarthy et al., 2020, and references therein). However, during the collisional stage, our thrust-mode models (for example model GC18, Fig. 10e–h) seem to reproduce some of the first-order features observed in the Pyrenees. We therefore suggest that the buoyancy push of subducted crust was insignificant during the Pyrenean orogeny.

Instead, convergence between Iberia and Europe induced shear-force-driven crustal wedging without deep (> 80 km) subduction of continental upper crust, as for example modelled by Jammes and Huismans (2012) and Grool et al. (2019).

### 4.3.2 Implications for the Alpine orogeny

In the Western Alps, a rifting phase prior to subduction lead to the formation of a ca. 300–400 km wide basin floored by exhumed mantle, and presumably only minor volumes of mature oceanic crust have been produced (Le Breton et al., 2021). Instead, ophiolites preserved in the Western Alps indicate serpentinization of the mantle exhumed in the basin (McCarthy et al., 2020, and references therein). The serpentinized material likely formed a relatively weak subduction interface (Zhao et al., 2020) and inhibited subduction of hydrous sediments, explaining the sparse arc magmatism in the European Alps (McCarthy et al., 2018, 2020; Yang et al., 2020). Field evidence of (ultra)high-pressure units (Chopin, 1984) in the Western Alps indicate either deep subduction of upper continental crust (Berger and Bousquet, 2008) at a close-to-lithostatic stress state, or significant deviation from the lithostatic stress state (Schenker et al., 2015) (Fig. 1b). Tomographic images (Zhao et al., 2015; Schmid et al., 2017) indicate that major volumes of upper continental crust have been involved in subduction. While the crustal wedge model may be applicable to the deformation during the post-collisional (≈30 Ma) stage of the Alpine orogeny (e.g. Erdős et al., 2019), it does not predict the exhumation
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Convergence and collisional stage of AC1 and AC18 (low and high shear resistance of serpentinite, complex density model, feldspar-dominated upper crust). White to dark red indicates the effective viscosity of the mantle material calculated by the algorithm, yellow to orange indicates the crustal matrix and the weak and strong inclusions of the overriding plate, light to dark blue indicates the corresponding crustal units of the subducting plate. The material parameters used for the upper crustal phases of the overriding and the subducting plates are identical. Green indicates the lower crust of both plates and salmon and brown indicates the sedimentary units. White lines are different isotherms and the arrows indicate the velocity field calculated by the numerical algorithm.</figcaption></figure>

tion of subduction-related (ultra)high-pressure continental and oceanic crustal rocks prior to collision. Instead, a subduction channel model has been proposed to explain deep subduction of continental upper crust and subsequent exhumation of (ultra)high-pressure units along the subduction interface (Gerya et al., 2002; Raimbourg et al., 2007; Butler et al., 2014). The subduction channel model, as proposed by Butler et al. (2014), has been criticized mainly for two reasons: (1) exhumation relies on significant removal, or erosion, of major crustal volumes, which is not in agreement with the sediment volume recorded in the Eocene–Oligocene basins (Malusà et al., 2015). (2) The exhuming units are strongly mixed (tectonic mélange) and significant volumes of lower crust are also exhumed, which is at odds with interpretations from seismic tomography showing no significant exhumation of lower crust (Schmid et al., 2017). In our model,

significant synconvergent exhumation of upper crust can occur by either diapirism or channel-flow and is enabled by spatially localized upper plate extension (Fig. 12) and not by significant erosion. Exhumation of the strong lower crust does not occur in our model, because it remains coupled to the subducting mantle lithosphere. Instead, weak zones in the upper crust connect and form a decoupling horizon within the upper crust above the Moho (see region of reduced viscosity in Fig. 7b). This observation is in agreement with interpretations from tomographic images. We focused here on the general evolution of different modes of orogenic wedges. We will present a detailed analysis of (i) the mechanisms causing local upper plate extension, (ii) the exhumation mechanisms of (ultra)high-pressure rocks, and (iii) the pressure and temperatures paths and associated exhumation velocities in a
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Wedge geometry at 170 Myr in model history. Blue to red is the density field calculated by the algorithm for (a) AS1 (LEOS) and (b) AC1 (CD model). White lines are several isotherms. (c) Topographic elevation of AS1 (red line) and AC1 (blue line).</figcaption></figure>

subsequent study, and we will not discuss these issues further here.

Our models are restricted to two dimensions and driven by far-field kinematic boundary conditions. Therefore, we can only capture first-order fundamental features of natural orogens. During model evolution, more and more crustal material is forced into subduction (Fig. 11d). In fact, in our models plate driving forces eventually reach again the magnitude necessary for subduction initiation after basin closure (see grey area in Fig. 11b). From this point on, the model probably becomes unrealistic, because in nature a new subduction zone may form at a different location and the active subduction zone may cease or slow down significantly. Plate reconstructions from the Western Alps indicate that European subduction below Adria was presumably slowed down significantly (“choked”) between ca. 35 and 25 Ma (Malusà et al., 2015; Schmid et al., 2017; see Fig. 13). Estimating subduction initiation at the Adriatic margin between ca. 90 and 85 Ma (Manzotti et al., 2014), there has been ca. 50–65 Myr of convergence before European subduction was “choked”. This duration approximately coincides with the time span (ca. 55–60 Myr of convergence) necessary to build up plate driving forces exceeding the magnitude for subduction initiation (at ca. 167 Myr; see 5 in Fig. 11b) in our model GC1. Hence, we suggest that the horizontal forces, required to drive the collision, became continuously larger during the evolving Alpine orogeny and ultimately exceeded a critical value leading to (1) significant slow-down of the European subduction and (2) plate boundary reorganization leading to subduction initiation of the remaining Adriatic oceanic lithosphere below Iberia (Fig. 13).

### 4.3.3 Challenges in applying wedge models to the Alpine orogeny

The “classical” crustal wedge models (Platt, 1986; Dahlen, 1990; Malavieille, 2010; Dal Zilio et al., 2020b) focus mainly on upper-crustal levels (see grey framed area in Fig. 1b and explanations in the introduction). The evolution of the Alpine orogeny, however, involves the entire lithosphere in wedging (Nicolas et al., 1990). Some studies applied the “classical” crustal wedge model to lithospheric-scale orogens (Platt, 1986; Vanderhaeghe et al., 2003; Beaumont et al., 2010). Vanderhaeghe et al. (2003) concluded that temperature dependence of material parameters indeed modulates wedge geometries, impacts on the transition between wedge and plateau behaviour, and controls the topographic evolution. In crustal wedge models, deformation is usually driven by an internal kinematic boundary condition at the base of the continental crust pulling the material towards a rigid backstop. Of course, subduction and collision in our models are also controlled by the kinematic boundary conditions. However, in our models boundary conditions are imposed far away from the evolving wedge, and not directly at the base of the continental crust, avoiding kinematic singularity points and allowing for continental subduction. This is likely important when the mantle lithosphere is involved in wedging (Beaumont et al., 2010). The first-order dynamics within the evolving wedge are controlled by the interaction of buoyancy and shear forces according to the shear resistance and density structure of the material. Our models, therefore, presumably provide a more realistic insight into orogen dynamics compared to orogenic wedge models considering crustal deformation only.
<figure><img src="image_12.png" /></figure>

**Figure 10.** Convergence and collisional stage of GC1 and GC18 (low and high shear resistance of serpentinite, CD model, quartz-dominated upper crust). White to dark red indicates the effective viscosity of the mantle material calculated by the algorithm, yellow to orange indicates the crustal matrix and the weak and strong inclusions of the overriding plate, light to dark blue indicates the corresponding crustal units of the subducting plate. The material parameters used for the upper crustal phases of the overriding and the subducting plates are identical. Green indicates the lower crust of both plates and salmon and brown indicates the sedimentary units. White lines are different isotherms and the arrows indicate the velocity field calculated by the numerical algorithm.

As mentioned in Sect. 4.1, the density structure impacts on the force balance and reduces the mean topographic elevation. This is important, because orogenic wedge models applied to the Western Alps have been criticized for overestimating the mean topographic elevation (Kissling, 1993; Dal Zilio et al., 2020a). Certainly, the average topography in our models is still higher compared to the average topography of the Alps. However, we suggest that employing more realistic density models is an important step to avoid exaggerated mean topographic elevation in orogenic wedge models. Furthermore, the topography in our models would likely decrease if we significantly reduce the convergence velocities at the mature stage of orogenic wedge formation, so that rollback of the subducted mantle lithosphere would become more important. The so-called rollback orogeny (Kissling

and Schlunegger, 2018; Dal Zilio et al., 2020a) applies to the latest Alpine evolution in the post-collisional stage (younger than ca. 30 Ma) but does not apply to the formation of the major (ultra)high-pressure regions (including oceanic, e.g. Zermatt-Saas, and Adriatic, e.g. Sesia and Dent Blanche, domains) of the Alps.

Absolute values for shear and buoyancy forces reported here strongly depend on the amount and strength of crustal material carried into the subduction zone. How much crustal material has been involved in Alpine subduction depends on the pre-orogenic crustal thickness (Mohn et al., 2014) and is still contentious (Butler, 2013; Schmid et al., 2017). The orogenic wedges modelled here are wider and deeper than the natural Alpine orogenic wedge. The size of the mechanical heterogeneities employed here is chosen based on the nu-
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Force evolution. (a) Ratio of buoyancy force to shear force, Ar<sub>F</sub>, (b) shear forces, (c) buoyancy forces, (d) average area of subducted material and (e) maximum depth of subducted crust. For a detailed explanation on the calculation of quantities see Appendix C. Numbering in all panels exemplified for model GC1: 1 = onset of subduction initiation, 2 = basin closure and onset of continental subduction, 3 = onset of continental crustal exhumation and first slab detachment, 4 = onset of the second exhumation event and 5 = onset of the second slab detachment.</figcaption></figure>

merical resolution (1 km × 1 km) and probably strongly im-
pacts on the size of the modelled wedge. At higher numerical
resolution, the size of the mechanical heterogeneities could
be reduced and the modelled rifted margins might be thinner,
leading to a more realistic pre-orogenic crustal thickness gen-
erated during the rifting period. In consequence, less material
would be involved in subduction, and the absolute magnitude
of forces would change. However, the general model evolu-
tion, i.e. the relative increase in buoyancy forces with con-
tinuous subduction, would be most likely unchanged. There-
fore, our models likely provide a representative model for
the relative evolution of shear and buoyancy forces building
orogenic wedges.

5 Conclusions

Our models show that upward-directed buoyancy forces, F<sub>B</sub>,
caused by subduction of continental crust, can be as high as
the horizontal driving shear forces, F<sub>D</sub>, induced by far-field
plate convergence. Therefore, such buoyancy forces should
be considered in models of continent collision and associated
orogenic wedge formation.

We investigated three parameters that control orogenic
wedge formation: (1) the shear resistance, or strength, of ser-
pentinites, controlling the strength of the subduction inter-
face; (2) the strength of the continental upper crust, control-
ling the maximal depth of crustal subduction; (3) the density
structure of the subducted material, which controls buoyancy
forces and significantly impacts on the mean topographic el-
evation, leading to more realistic topography evolution.

These three parameters control the evolution of the ra-
tio F<sub>B</sub>/F<sub>D</sub> = Ar<sub>F</sub>. The value of Ar<sub>F</sub> controls the mode of
<figure><img src="image_14.png" /></figure>

**Figure 12.** Orogeny dynamics observed in the models presented here. Light blue and orange indicates the upper crustal phases, dark blue and merlot red indicates the lower crustal phases of the subducting and upper plates, respectively. Dark grey and magenta represent the sediments and the serpentinite unit, respectively. Arrows indicate the velocity field calculated by the numerical algorithm.

<figure><img src="image_15.png" /></figure>

**Figure 13.** Simplified plate reconstruction of Europe, Iberia and Adria at (a) 35 Ma and (b) 25 Ma. Europe and Iberia are fixed, Adria is moving northward. Thick black and grey solid lines indicate active subduction zones, dashed lines indicate inactive subduction zones. Sketches modified after Schmid et al. (2017) and Malusà et al. (2015).

orogenic wedge formation: $Ar_F \approx 0.5$ causes thrust-sheet-dominated wedges, similar to crustal wedge models without buoyancy, with thrusts propagating towards the foreland; $Ar_F \approx 0.75$ causes minor wedge formation due to significant relamination of subducted crust below the upper plate; $Ar_F \approx 1$ causes buoyancy flow, or diapir, dominated wedges involving exhumation of crustal material from large depth (> 80 km).

The spatial and temporal variation of the three parameters mentioned above, and the associated variation of $Ar_F$, may explain the variation of deformation style observed between and within natural collisional orogens. We suggest that during the formation of the Pyrenees $Ar_F \lesssim 0.5$ due to the absence of high-grade metamorphic rocks, whereas for the Alps $Ar_F \approx 1$ during exhumation of high-grade rocks and $Ar_F \lesssim 0.5$ during the post-collisional stage. In the models, the
increase in $F_B$ during orogenic wedge growth causes an increase in $F_D$ required to drive wedge growth and subduction. In nature, the requirement of larger driving forces may cause a slow down, or “choking”, of the associated subduction and may even cause horizontally forced subduction initiation in other, nearby regions. Therefore, quantifying buoyancy and shear forces during orogenic wedge formation may prove useful to unravel changes in relative plate motion and subduction initiation during the Alpine orogeny, and during other collisional orogenies worldwide.
## Appendix A: Algorithm description

We employ the extended Boussinesq approximation (e.g. Candioti et al., 2020) for buoyancy-driven flow. The applied numerical algorithm solves the continuity and momentum equations defined as

$$ \frac{\partial v_i}{\partial x_i} = 0, \qquad (A1) $$

$$ \frac{\partial \sigma_{ij}}{\partial x_j} = -\rho g_i, \qquad (A2) $$

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 the gravitational acceleration. Effective density can be (1) computed as a linearized equation of state (LEOS) like

$$ \rho(P, T) = \rho_0 (1 - \alpha \Delta T) (1 + \beta \Delta P), \qquad (A3) $$

where P (negative mean stress) is pressure, T is temperature, ρ₀ is the material density at reference temperature (T₀) and pressure (P₀), and α and β are material parameters accounting for density changes due to thermal expansion and isothermal compression, respectively, ΔT = T - T₀, and ΔP = P - P₀. Alternatively, effective density can be (2) pre-computed using the Gibbs free energy minimization software package Perple_X (Connolly, 2005) for given bulk rock compositions. The density field predicted by these phase equilibrium models is stored in look-up tables and density values are read in by the applied algorithm during simulation run-time according to local pressure and temperature conditions at each grid cell centre.

We employ a backward-Euler scheme (e.g. Schmalholz et al., 2001) to define the viscoelastic stress tensor components as

$$ \sigma_{ij} = -P\delta_{ij} + 2 \eta^{\text{eff}} \dot{\epsilon}_{ij}^{\text{eff}}, \qquad (A4) $$

$$ \delta_{ij} = \begin{cases} 0, & i \neq j, \\ 1, & i = j, \end{cases} \qquad (A5) $$

where δij is the Kronecker-Delta function, ηeff is the effective viscosity, and εijeff denotes the effective deviatoric strain rate tensor components,

$$ \dot{\epsilon}_{ij}^{\text{eff}} = \left( \dot{\epsilon}_{ij} + \frac{\tau_{ij}^{o}}{2G\Delta t} \right), \qquad (A6) $$

where G is shear modulus, Δt is the time step, τij denotes the deviatoric stress tensor components of the preceding time step. We consider Maxwell materials and additively decompose the total deviatoric strain rate tensor εij into contributions from viscous (dislocation, diffusion and Peierls creep), elastic and plastic deformation like

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{\text{ela}} + \dot{\epsilon}_{ij}^{\text{pla}} + \dot{\epsilon}_{ij}^{\text{dis}} + \dot{\epsilon}_{ij}^{\text{dif}} + \dot{\epsilon}_{ij}^{\text{pei}}. \qquad (A7) $$

Furthermore, we perform an iteration cycle locally on each grid cell until Eq. (A7) is satisfied (e.g. Popov and Sobolev, 2008). The viscosity for the dislocation and Peierls creep flow laws is a function of the second invariant of the respective strain rate components, and the dislocation creep viscosity takes the following form:

$$ \eta^{\text{dis}} = \frac{2^{\frac{1-n}{n}}}{3^{\frac{1+n}{2n}}} \zeta A^{-\frac{1}{n}} (\dot{\epsilon}_{\text{II}}^{\text{dis}})^{\frac{1}{n}-1} \exp\left(\frac{Q+PV}{nRT}\right) (f_{\text{H}_2\text{O}})^{-\frac{r}{n}}, \qquad (A8) $$

where the ratio in front of the prefactor ζ results from tensor conversion of the experimentally derived flow law (e.g. Schmalholz and Fletcher, 2011). A, n, Q, V, fH₂O and r are material parameters. The diffusion creep viscosity is calculated as

$$ \eta^{\text{dif}} = \frac{1}{3} A^{-1} d^m \exp\left(\frac{Q+PV}{RT}\right) (f_{\text{H}_2\text{O}})^{-r}, \qquad (A9) $$

where d is grain size and m is a grain size exponent. Effective Peierls viscosity is calculated based on the experimentally derived flow law by Goetze and Evans (1979) expressed in the regularized form of Kameyama et al. (1999) as

$$ \eta^{\text{pei}} = \frac{2^{\frac{1-s}{s}}}{3^{\frac{1+s}{2s}}} \hat{A} (\dot{\epsilon}_{\text{II}}^{\text{pei}})^{\frac{1}{s}-1}, \qquad (A10) $$

where s is an effective stress exponent:

$$ s = 2 \gamma \frac{Q}{RT} (1 - \gamma). \qquad (A11) $$

$\hat{A}$ in Eq. (A10) is

$$ \hat{A} = \left[ A_P \exp\left(-\frac{Q(1-\gamma)^2}{RT}\right) \right]^{-\frac{1}{s}} \gamma \sigma_P, \qquad (A12) $$

where A_P, γ and σ_P are flow law parameters. Brittle-plastic material failure is controlled by the Drucker-Prager yield function

$$ F = \tau_{\text{II}} - P \sin\phi - C \cos\phi, \qquad (A13) $$

which depends on the internal friction angle, φ, and the cohesion, C. In case of failure (F ≥ 0), the plastic viscosity at the yield stress is calculated as

$$ \eta^{\text{pla}} = \frac{P \sin\phi + C \cos\phi}{2\dot{\epsilon}_{\text{II}}^{\text{eff}}}, \qquad (A14) $$

and εIIeff = √(1/2 εijeff εijeff). In Eq. (A4), the effective viscosity is either the quasi-harmonic average of the visco-elastic contributions

$$ \eta^{\text{eff}} = \begin{cases} \left( \frac{1}{G\Delta t} + \frac{1}{\eta^{\text{dis}}} + \frac{1}{\eta^{\text{dif}}} + \frac{1}{\eta^{\text{pei}}} \right)^{-1} & , F < 0, \\ \eta^{\text{pla}} & , F \ge 0, \end{cases} \qquad (A15) $$
<table><caption>**Table A1.** Physical parameters used in the numerical simulations presented here.</caption><thead><tr><th>Model unit</th><th>Rheology (reference)</th><th>k [W m<sup>-1</sup> K<sup>-1</sup>]</th><th>H<sub>R</sub> [W m<sup>-3</sup>]</th><th>C [Pa]</th><th>φ [°]</th></tr></thead><tbody><tr><td>Crustal matrix 1<sup>a,b</sup></td><td>Wet anorthite (Rybacki and Dresen, 2004)</td><td>2.25</td><td>1.0200 × 10<sup>-6</sup></td><td>1 × 10<sup>7</sup></td><td>30</td></tr><tr><td>Crustal matrix 2<sup>a,b</sup></td><td>Westerly granite (Hansen and Carter, 1983)</td><td>2.25</td><td>1.0200 × 10<sup>-6</sup></td><td>1 × 10<sup>7</sup></td><td>30</td></tr><tr><td>Weak inclusion<sup>a,b</sup></td><td>Wet quartzite (Ranalli, 1995)</td><td>2.25</td><td>1.0200 × 10<sup>-6</sup></td><td>1 × 10<sup>6</sup></td><td>5</td></tr><tr><td>Strong inclusion<sup>a,b</sup></td><td>Maryland diabase (Mackwell et al., 1998)</td><td>2.25</td><td>1.0200 × 10<sup>-6</sup></td><td>1 × 10<sup>7</sup></td><td>30</td></tr><tr><td>Calcite<sup>a,b</sup></td><td>Calcite (Schmid et al., 1977)</td><td>2.37</td><td>0.5600 × 10<sup>-6</sup></td><td>1 × 10<sup>7</sup></td><td>30</td></tr><tr><td>Mica<sup>a,b</sup></td><td>Mica (Kronenberg et al., 1990)</td><td>2.55</td><td>2.9000 × 10<sup>-6</sup></td><td>1 × 10<sup>7</sup></td><td>15</td></tr><tr><td>Lower crust<sup>a,c</sup></td><td>Wet anorthite (Rybacki and Dresen, 2004)</td><td>2.25</td><td>0.2600 × 10<sup>-6</sup></td><td>1 × 10<sup>7</sup></td><td>30</td></tr><tr><td>Strong mantle<sup>a,d</sup></td><td>Dry olivine (Hirth and Kohlstedt, 2003)</td><td>2.75</td><td>2.1139 × 10<sup>-8</sup></td><td>1 × 10<sup>7</sup></td><td>30</td></tr><tr><td>Weak mantle<sup>a,d</sup></td><td>Wet olivine (Hirth and Kohlstedt, 2003)</td><td>2.75</td><td>2.1139 × 10<sup>-8</sup></td><td>1 × 10<sup>7</sup></td><td>30</td></tr><tr><td>Serpentinite<sup>a,e</sup></td><td>Antigorite (Hilairet et al., 2007)</td><td>2.75</td><td>2.1139 × 10<sup>-8</sup></td><td>1 × 10<sup>7</sup></td><td>25</td></tr><tr><td>Dislocation creep</td><td>A [Pa<sup>-n-r</sup> s<sup>-1</sup>]</td><td>ζ [ ]</td><td>n [ ]</td><td>Q [J mol<sup>-1</sup>]</td><td>V [m<sup>3</sup> mol<sup>-1</sup>]</td><td>r [ ]</td></tr><tr><td>Crustal matrix 1</td><td>3.9811 × 10<sup>-16</sup></td><td>0.3</td><td>3.0</td><td>356 × 10<sup>3</sup></td><td>0.00 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Crustal matrix 2</td><td>3.1623 × 10<sup>-26</sup></td><td>1.0</td><td>3.3</td><td>186.5 × 10<sup>3</sup></td><td>0.00 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Weak inclusion</td><td>5.0717 × 10<sup>-18</sup></td><td>1.0</td><td>2.3</td><td>154 × 10<sup>3</sup></td><td>0.00 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Strong inclusion</td><td>5.0477 × 10<sup>-28</sup></td><td>1.0</td><td>4.7</td><td>485 × 10<sup>3</sup></td><td>0.00 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Calcite</td><td>1.5849 × 10<sup>-25</sup></td><td>1.0</td><td>4.7</td><td>297 × 10<sup>3</sup></td><td>0.00 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Mica</td><td>1.0000 × 10<sup>-138</sup></td><td>1.0</td><td>18.0</td><td>51.0 × 10<sup>3</sup></td><td>0.00 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Lower crust</td><td>3.9811 × 10<sup>-16</sup></td><td>1.0</td><td>3.0</td><td>356 × 10<sup>3</sup></td><td>0.00 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Strong mantle</td><td>1.1000 × 10<sup>-16</sup></td><td>1.0</td><td>3.5</td><td>530 × 10<sup>3</sup></td><td>14.0 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Weak mantle<sup>f</sup></td><td>5.6786 × 10<sup>-27</sup></td><td>1.0</td><td>3.5</td><td>480 × 10<sup>3</sup></td><td>11.0 × 10<sup>-6</sup></td><td>1.2</td></tr><tr><td>Serpentinite</td><td>4.4738 × 10<sup>-38</sup></td><td>see Table 1</td><td>3.8</td><td>8.90 × 10<sup>3</sup></td><td>3.20 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Diffusion creep<sup>g</sup></td><td>A [Pa<sup>-n-r</sup> m<sup>m</sup> s<sup>-1</sup>]</td><td>m [ ]</td><td>n [ ]</td><td>Q [J mol<sup>-1</sup>]</td><td>V [m<sup>3</sup> mol<sup>-1</sup>]</td><td>r [ ]</td></tr><tr><td>Strong mantle</td><td>1.5000 × 10<sup>-15</sup></td><td>3.0</td><td>1.0</td><td>370 × 10<sup>3</sup></td><td>7.5 × 10<sup>-6</sup></td><td>0.0</td></tr><tr><td>Weak mantle<sup>f</sup></td><td>2.5000 × 10<sup>-23</sup></td><td>3.0</td><td>1.0</td><td>375 × 10<sup>3</sup></td><td>9.0 × 10<sup>-6</sup></td><td>1.0</td></tr><tr><td>Peierls creep</td><td>A<sub>P</sub> [s<sup>-1</sup>]</td><td>Q [J mol<sup>-1</sup>]</td><td>V [m<sup>3</sup> mol<sup>-1</sup>]</td><td>σ<sub>P</sub> [Pa]</td><td>γ [ ]</td><td></td></tr><tr><td>Mantle<sup>h</sup></td><td>5.7000 × 10<sup>11</sup></td><td>540 × 10<sup>3</sup></td><td>0.0 × 10<sup>-6</sup></td><td>8.5 × 10<sup>9</sup></td><td>0.1</td><td></td></tr></tbody></table>

Constant parameters: a a heat capacity $c_p = 1050$ [J kg$^{-1}$ K$^{-1}$] is employed in Eq. (A20) for all phases. b,c,d A constant shear modulus $G = 2 \times 10^{10}$ [Pa] is used in Eq. (A4). e A constant shear modulus $G = 1.81 \times 10^{10}$ [Pa] is used in Eq. (A4). b $\rho_0 = 2800$ [kg m$^{-3}$] and e $\rho_0 = 2900$ [kg m$^{-3}$] and b,c $\alpha = 3.5 \times 10^{-5}$ [K$^{-1}$] and e $\alpha = 4.7 \times 10^{-5}$ [K$^{-1}$] and b,c,e $\beta = 1 \times 10^{-11}$ [Pa$^{-1}$] are used for density calculations using a simplified equation of state according to Eq. (A4). f A water fugacity $f_{H_2O} = 1.0 \times 10^9$ [Pa] is used in equation Eq. (A8).

For all other phases $f_{H_2O} = 0.0$ [Pa]. g A constant grain size $d = 1 \times 10^{-3}$ [m] is used in Eq. (A9). h Reference: Goetze and Evans (1979) regularized by Kameyama et al. (1999). These parameters are used for both strong and weak mantle rheology.

or is equal to the viscosity $\eta^{\text{pla}}$ calculated according to Eq. (A14). Rigid body rotation is computed analytically:

$$ \tau_{ij} = \mathbf{R}^{\text{T}} \tau_{ij} \mathbf{R}, \qquad (\text{A16}) $$

$$ \mathbf{R} = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix}, \qquad (\text{A17}) $$

$$ \theta = \Delta t \omega_{ij}, \qquad (\text{A18}) $$

$$ \omega_{ij} = \frac{1}{2} \left( \frac{\partial v_j}{\partial x_i} - \frac{\partial v_i}{\partial x_j} \right), \qquad (\text{A19}) $$

where $\mathbf{T}$ denotes matrix transposition, $\mathbf{R}$ is the rotation matrix, $\theta$ is the rotation angle and $\omega_{ij}$ denotes vorticity tensor components. Heat transfer is included into the model by expressing the energy balance equation with respect to temperature as

$$ \rho c_p \frac{DT}{Dt} = \frac{\partial}{\partial x_i} \left( k \frac{\partial T}{\partial x_i} \right) + H_A + H_D + H_R, \qquad (\text{A20}) $$

where $c_p$ is the specific heat capacity, $D/Dt$ is the material derivative, $k$ is thermal conductivity, $H_A = T\alpha v_z g\rho$ includes

contributions from adiabatic processes assuming lithostatic pressure conditions, $H_D = \tau_{ij}^2/(2\eta^{\text{eff}})$ includes contributions from dissipative processes and $H_R$ includes heat production from radiogenic elements.

## Appendix B: Emplacement of elliptical heterogeneous inclusions

The semi-major axis is $a = 30$ km and the semi-minor axis is $b = 2.5$ km for all elliptical inclusions. The inclusions are emplaced between $-400$ km $< x < 400$ km at two different vertical levels. The $x$ coordinate of the first elliptical inclusion's centre $x_n^{\text{C,L}}$ (superscript C = crust, L = lithosphere) at each vertical level is calculated as

$$ x_n^{\text{C,L}} = (R^\text{E} - R^\text{S}) A + R^\text{S}, \qquad (\text{B1}) $$

$$ R^\text{S} = (d - c) A - d, \qquad (\text{B2}) $$

$$ R^\text{E} = (f - h) A - f, \qquad (\text{B3}) $$
where $R^S$ is a random starting $x$ coordinate, $R^E$ is a random ending $x$ coordinate, $d = 400$ km, $c = 350$ km, $0 < A < 1$ is a random amplitude, $f = 300$ km and $h = 250$ km. The starting coordinate of the next horizontally emplaced inclusion $x_{n+1}^{C,L}$ is then calculated as

$$x_{n+1}^{C,L} = x_n^{C,L} + R^{dx}, \qquad (B4)$$

$$R^{dx} = a (1.75 + A), \qquad (B5)$$

where $R^{dx}$ is a random spacing ensuring that the ellipses do not overlap each other. The $z$ coordinate of the elliptical inclusion's centre $z_m^C$ at the vertical level $m$ is calculated as

$$z_1^C = \left(\frac{1}{3}z_{\text{Moho}} - 3b\right) A - \frac{1}{3}z_{\text{Moho}} + b, \qquad (B6)$$

$$z_2^C = \left(\frac{1}{3}z_{\text{Moho}} - 2b\right) A - \frac{2}{3}z_{\text{Moho}} + b, \qquad (B7)$$

where $z_{\text{Moho}} = 33$ km is the initial depth of the Moho. Similarly, the central coordinates of the elliptical inclusions in the mantle are calculated using the following parameters for the first level:

$$R^S = (d-c) A - d, \qquad (B8)$$

$$R^E = (f-h) A - f, \qquad (B9)$$

$$R^{dx} = a (4 + A), \qquad (B10)$$

$$z_1^L = \frac{1}{3} \left( z_{\text{Lith}} - z_{\text{Moho}} - 9b \right) A - \frac{1}{3} \left( z_{\text{Lith}} - z_{\text{Moho}} \right) - z_{\text{Moho}} + b, \qquad (B11)$$

and for the second level:

$$R^S = (d-c) A - d - 5a, \qquad (B12)$$

$$R^E = (f-h) A - f, \qquad (B13)$$

$$R^{dx} = a (4 + A), \qquad (B14)$$

$$z_2^L = \frac{1}{3} \left( z_{\text{Lith}} - z_{\text{Moho}} \right) A - z_{\text{Lith}}, \qquad (B15)$$

where $z_{\text{Lith}} = 60$ km. Finally, all markers within the circumference of the elliptical inclusion are assigned with a random phase, i.e. either mechanically weak or strong material, according to the condition

$$\frac{(x_M - x_n^{C,L})^2}{a^2} + \frac{(z_M - z_m^{C,L})^2}{b^2} < 1, \qquad (B16)$$

where $x_M$ and $z_M$ are the horizontal and vertical coordinate of the marker, respectively. All random numbers used here are seeded at a value of 197 using the C function `srand` (Intel compiler version 13.1.3). Choosing the above-mentioned values yields an increased number of weak elliptical inclusions in the centre of the domain. This yields localization of deformation without an additional perturbation of the marker field in the centre of the domain.

## Appendix C: Buoyancy and driving forces

In this study we use $F_D = 2 \times \bar{\tau}_{\text{II}}^{\text{avg}}$ as a representative value for the horizontal driving force per unit length causing the far-field convergence. This horizontal driving force calculation is correct if the smallest principal stress is vertical and corresponds to the lithostatic pressure and if the maximal principal stress is horizontal (e.g. Schmalholz et al., 2014, 2019), which is the case around the lateral model boundaries where $F_D$ is calculated. To calculate $\bar{\tau}_{\text{II}}^{\text{avg}}$, first the vertical integral of $\tau_{\text{II}}(x, y)$ is calculated at the different horizontal positions of the numerical grid:

$$\bar{\tau}_{\text{II}}(x) = \int_{\text{Sb}}^{\text{St}(x)} \tau_{\text{II}}(x, z) \, dz. \qquad (C1)$$

The values of $\bar{\tau}_{\text{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 $\bar{\tau}_{\text{II}}^{\text{avg}}$. $\text{St}(x)$ is the height of the topography at each horizontal grid point and $\text{Sb} = -660$ km is the bottom of the domain. The reader is referred to Candioti et al. (2020) for further detail.

The buoyancy force per unit length is calculated with the difference between densities of all subducted material except for mantle (i.e. upper and lower crust, serpentinite, and sediments) and a representative mantle density of 3350 kg m⁻³. This mantle density is representative for mantle rocks down to depths of ≈ 150 km (Candioti et al., 2020), which is, in our models, a representative maximal depth of subduction of non-mantle material. The buoyancy force is calculated by

$$F_B = \int_{\Omega} \Delta \rho g \, d\Omega, \qquad (C2)$$

where $\Delta\rho$ is the above-mentioned density difference, which is integrated over the area $\Omega$ occupied by non-mantle material subducted below $z = -40$ km. To obtain a reasonable value for the upward-directed buoyancy force of subducted material, which is not isostatically balanced by high topography, we subtract the force contributing to the build-up of topography. As an approximation for this topographic contribution, the density of material lifted above $z = 1.5$ km is integrated over its area and then subtracted from $F_B$. As a measure for the relative importance of buoyancy and shear forces in orogen dynamics, we here define $Ar_F = F_B/F_D$ which is a dimensionless number comparable to the Argand number (England and McKenzie, 1982).
*Data availability.* The data presented in this study are available on request from Lorenzo G. Candioti.

*Video supplement.* We provide videos showing the entire evolution of the numerical simulation REF (https://doi.org/10.5446/50527, Candioti, 2020a) and the subduction and collision stage of the numerical simulation GC1 (https://doi.org/10.5446/50528, Candioti, 2020b).

*Author contributions.* LGC configured and performed the numerical simulations, interpreted the numerical results, generated the figures and wrote the paper. SMS designed the scientific and numerical study, helped in interpreting the results and designing the figures, contributed to writing the paper, and acquired the funding for this study. TD developed the applied numerical algorithm and helped in configuring the model and the interpretation of the results. EM performed the phase equilibria calculations and helped interpreting the results.

*Competing interests.* The authors declare that they have no conflict of interest.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

*Special issue statement.* This article is part of the special issue "New insights into the tectonic evolution of the Alps and the adjacent orogens". It is not associated with a conference.

*Acknowledgements.* The authors gratefully thank Jianfeng Yang and Jonas Ruh for very constructive comments during the review process.

*Financial support.* This research has been supported by the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (grant no. 200020 163169).

*Review statement.* This paper was edited by Mirijam Vrabec and reviewed by Jonas B. Ruh and Jianfeng Yang.

## References

* Austin, N. J. and Evans, B.: Paleowattmeters: A scaling relation for dynamically recrystallized grain size, Geology, 35, 343–346, 2007.

* Austrheim, H.: Eclogitization of lower crustal granulites by fluid migration through shear zones, Earth Planet. Sc. Lett., 81, 221–232, 1987.

* Auzemery, A., Willingshofer, E., Yamato, P., Duretz, T., and Sokoutis, D.: Strain localization mechanisms for subduction initiation at passive margins, Global Planet. Change, 195, 103323, https://doi.org/10.1016/j.gloplacha.2020.103323, 2020.

* Barnhoorn, A., Drury, M. R., and van Roermund, H. L.: Evidence for low viscosity garnet-rich layers in the upper mantle, Earth Planet. Sc. Lett., 289, 54–67, 2010.

* Bauville, A. and Schmalholz, S. M.: Transition from thin-to thick-skinned tectonics and consequences for nappe formation: Numerical simulations and applications to the Helvetic nappe system, Switzerland, Tectonophysics, 665, 101–117, 2015.

* Beaumont, C., Ellis, S., Hamilton, J., and Fullsack, P.: Mechanical model for subduction-collision tectonics of Alpine-type compressional orogens, Geology, 24, 675–678, 1996.

* Beaumont, C., Jamieson, R., and Nguyen, M.: Models of large, hot orogens containing a collage of reworked and accreted terranes, Can. J. Earth Sci., 47, 485–515, 2010.

* Behr, W. M. and Becker, T. W.: Sediment control on subduction plate speeds, Earth Planet. Sc. Lett., 502, 166–173, 2018.

* Berger, A. and Bousquet, R.: Subduction-related metamorphism in the Alps: review of isotopic ages based on petrology and their geodynamic consequences, Geol. Soc. Lond. Spec. Publ., 298, 117–144, 2008.

* Bessat, A., Duretz, T., Hetényi, G., Pilet, S., and Schmalholz, S. M.: Stress and deformation mechanisms at a subduction zone: insights from 2D thermo-mechanical numerical modelling, Geophys. J. Int., 221, 1605–1625, https://doi.org/10.1093/gji/ggaa092, 2020.

* Borderie, S., Graveleau, F., Witt, C., and Vendeville, B. C.: Impact of an interbedded viscous décollement on the structural and kinematic coupling in fold-and-thrust belts: Insights from analogue modeling, Tectonophysics, 722, 118–137, 2018.

* Bürgmann, R. and Dresen, G.: Rheology of the lower crust and upper mantle: Evidence from rock mechanics, geodesy, and field observations, Annu. Rev. Earth Pl. Sc., 36, 531–567, https://doi.org/10.1146/annurev.earth.36.031207.124326, 2008.

* Burov, E. and Watts, A.: The long-term strength of continental lithosphere: "jelly sandwich" or "crème brûlée"?, GSA Today, 16, 4, https://doi.org/10.1130/1052-5173(2006)016<4:tltSOc>2.0.cO;2, 2006.

* Burov, E., François, T., Agard, P., Le Pourhiet, L., Meyer, B., Tirel, C., Lebedev, S., Yamato, P., and Brun, J.-P.: Rheological and geodynamic controls on the mechanisms of subduction and HP/UHP exhumation of crustal rocks during continental collision: Insights from numerical models, Tectonophysics, 631, 212–250, 2014.

* Butler, J. P., Beaumont, C., and Jamieson, R. A.: The Alps 2: Controls on crustal subduction and (ultra) high-pressure rock exhumation in Alpine-type orogens, J. Geophys. Res.-Sol. Ea., 119, 5987–6022, 2014.

* Butler, R. W.: Area balancing as a test of models for the deep structure of mountain belts, with specific reference to the Alps, J. Struct. Geol., 52, 2–16, 2013.

* Byerlee, J.: Friction of rocks, in: Rock friction and earthquake prediction, 615–626, Springer, 1978.

* Candioti, L. G.: Evolution of numerical simulation REF, TIB, https://doi.org/10.5446/50527, 2020a.

* Candioti, L. G.: Evolution of numerical simulation GC1, TIB, https://doi.org/10.5446/50528, 2020b.
Candioti, L. G., Schmalholz, S. M., and Duretz, T.: Impact of upper mantle convection on lithosphere hyperextension and subsequent horizontally forced subduction initiation, Solid Earth, 11, 2327–2357, https://doi.org/10.5194/se-11-2327-2020, 2020.

Chapple, W. M.: Mechanics of thin-skinned fold-and-thrust belts, Geol. Soc. Am. Bull., 89, 1189–1198, 1978.

Chenin, P., Manatschal, G., Picazo, S., Müntener, O., Karner, G., Johnson, C., and Ulrich, M.: Influence of the architecture of magma-poor hyperextended rifted margins on orogens produced by the closure of narrow versus wide oceans, Geosphere, 13, 559–576, 2017.

Chenin, P., Picazo, S., Jammes, S., Manatschal, G., Müntener, O., and Karner, G.: Potential role of lithospheric mantle composition in the Wilson cycle: a North Atlantic perspective, Geol. Soc. Lond. Spec. Publ., 470, 157–172, 2019.

Chernak, L. J. and Hirth, G.: Deformation of antigorite serpentinite at high temperature and pressure, Earth Planet. Sc. Lett., 296, 23–33, 2010.

Chopin, C.: Coesite and pure pyrope in high-grade blueschists of the Western Alps: a first record and some consequences, Contrib. Mineral. Petr., 86, 107–118, 1984.

Connolly, J. A.: Computation of phase equilibria by linear programming: a tool for geodynamic modeling and its application to subduction zone decarbonation, Earth Planet. Sc. Lett., 236, 524–541, 2005.

Crameri, F.: Geodynamic diagnostics, scientific visualisation and StagLab 3.0, Geosci. Model Dev., 11, 2541–2562, https://doi.org/10.5194/gmd-11-2541-2018, 2018.

Crameri, F., Magni, V., Domeier, M., Shephard, G. E., Chotalia, K., Cooper, G., Eakin, C. M., Grima, A. G., Gürer, D., Király, Á., Mulyukova, E., Peters, K., Robert, B., and Thielmann, M.: A transdisciplinary and community-driven database to unravel subduction zone initiation, Nat. Commun., 11, 1–14, 2020.

Currie, C. A., Beaumont, C., and Huismans, R. S.: The fate of subducted sediments: A case for backarc intrusion and underplating, Geology, 35, 1111–1114, 2007.

Dahlen, F.: Critical taper model of fold-and-thrust belts and accretionary wedges, Annu. Rev. Earth Planet. Sc., 18, 55–99, 1990.

Dahlen, F., Suppe, J., and Davis, D.: Mechanics of fold-and-thrust belts and accretionary wedges: Cohesive Coulomb theory, J. Geophys. Res.-Sol. Ea., 89, 10087–10101, 1984.

Dal Zilio, L., Kissling, E., Gerya, T., and van Dinther, Y.: Slab Roll-back Orogeny model: A test of concept, Geophys. Res. Lett., 47, e2020GL089917, https://doi.org/10.1029/2020GL089917, 2020a.

Dal Zilio, L., Ruh, J., and Avouac, J.-P.: Structural evolution of orogenic wedges: interplay between erosion and weak décollements, Tectonics, 39, e2020TC006210, https://doi.org/10.1029/2020TC006210, 2020b.

Dannberg, J., Eilon, Z., Faul, U., Gassmöller, R., Moulik, P., and Myhill, R.: The importance of grain size to mantle dynamics and seismological observations, Geochem. Geophy. Geosy., 18, 3034–3061, 2017.

David, E. C., Brantut, N., Hansen, L. N., and Mitchell, T. M.: Absence of stress-induced anisotropy during brittle deformation in antigorite serpentinite, J. Geophys. Res.-Sol. Ea., 123, 10616–10644, 2018.

Duretz, T. and Gerya, T.: Slab detachment during continental collision: Influence of crustal rheology and interaction with lithospheric delamination, Tectonophysics, 602, 124–140, 2013.

Duretz, T., Schmalholz, S., and Gerya, T.: Dynamics of slab detachment, Geochem. Geophy. Geosy., 13, Q03020, https://doi.org/10.1029/2011GC004024, 2012.

Duretz, T., May, D. A., and Yamato, P.: A free surface capturing discretization for the staggered grid finite difference scheme, Geophys. J. Int., 204, 1518–1530, 2016a.

Duretz, T., Petri, B., Mohn, G., Schmalholz, S., Schenker, F., and Müntener, O.: The importance of structural softening for the evolution and architecture of passive margins, Sci. Rep.-UK, 6, 38704, https://doi.org/10.1038/srep38704, 2016b.

England, P. and McKenzie, D.: A thin viscous sheet model for continental deformation, Geophys. J. Int., 70, 295–321, 1982.

Erdős, Z., Huismans, R. S., van der Beek, P., and Thieulot, C.: Extensional inheritance and surface processes as controlling factors of mountain belt structure, J. Geophys. Res.-Sol. Ea., 119, 9042–9061, 2014.

Erdős, Z., Huismans, R. S., and van der Beek, P.: Control of increased sedimentation on orogenic fold-and-thrust belt structure – insights into the evolution of the Western Alps, Solid Earth, 10, 391–404, https://doi.org/10.5194/se-10-391-2019, 2019.

Forsyth, D. and Uyeda, S.: On the relative importance of the driving forces of plate motion, Geophys. J. Int., 43, 163–200, 1975.

Gerya, T.: Introduction to numerical geodynamic modelling, Cambridge University Press, 2019.

Gerya, T. V. and Yuen, D. A.: Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties, Phys. Earth Planet. In., 140, 293–318, 2003.

Gerya, T. V., Stöckhert, B., and Perchuk, A. L.: Exhumation of high-pressure metamorphic rocks in a subduction channel: A numerical simulation, Tectonics, 21, 6–1, 2002.

Gerya, T. V., Perchuk, L. L., Maresch, W. V., and Willner, A. P.: Inherent gravitational instability of hot continental crust: Implications for doming and diapirism in granulite facies terrains, SPECIAL PAPERS-GEOLOGICAL SOCIETY OF AMERICA, 97–116, 2004.

Goetze, C. and Evans, B.: Stress and temperature in the bending lithosphere as constrained by experimental rock mechanics, Geophys. J. Int., 59, 463–478, 1979.

Graveleau, F., Malavieille, J., and Dominguez, S.: Experimental modelling of orogenic wedges: A review, Tectonophysics, 538, 1–66, 2012.

Grool, A. R., Huismans, R. S., and Ford, M.: Salt décollement and rift inheritance controls on crustal deformation in orogens, Terra Nova, 31, 562–568, 2019.

Guillot, S., Schwartz, S., Reynard, B., Agard, P., and Prigent, C.: Tectonic significance of serpentinites, Tectonophysics, 646, 1–19, 1995.

Gutscher, M.-A., Kukowski, N., Malavieille, J., and Lallemand, S.: Episodic imbricate thrusting and underthrusting: Analog experiments and mechanical analysis applied to the Alaskan accretionary wedge, J. Geophys. Res.-Sol. Ea., 103, 10161–10176, 1998.

Hacker, B. R., Peacock, S. M., Abers, G. A., and Holloway, S. D.: Subduction factory 2. Are intermediate-depth earthquakes in subducting slabs linked to metamorphic dehy-
dration reactions?, J. Geophys. Res.-Sol. Ea., 108, 2030,
https://doi.org/10.1029/2001JB001129, 2003.

Handy, M. R., Schmid, S. M., Bousquet, R., Kissling, E., and Bernoulli, D.: Reconciling plate-tectonic reconstructions of Alpine Tethys with the geological-geophysical record of spreading and subduction in the Alps, Earth-Sci. Rev., 102, 121–158, 2010.

Hansen, F. and Carter, N.: Semibrittle creep of dry and wet Westerly granite at 1000 MPa, ARMA-83-0429, in: 24th U.S. Symposium on Rock Mechanics (USRMS), College Station, Texas, 1983.

Hansen, L. N., David, E. C., Brantut, N., and Wallis, D.: In-sight into the microphysics of antigorite deformation from spher-ical nanoindentation, Philos. T. R. Soc. A, 378, 20190197, https://doi.org/10.1098/rsta.2019.0197, 2020.

Henry, P., Azambre, B., Montigny, R., Rossy, M., and Steven-
son, R.: Late mantle evolution of the Pyrenean sub-continental
lithospheric mantle in the light of new 40Ar-39Ar and Sm-Nd
ages on pyroxenites and peridotites (Pyrenees, France), Tectono-
physics, 296, 103–123, 1998.

Hess, H. H.: Serpentines, orogeny, and epeirogeny, Geol. Soc. Am. Spec. Paper, 62, 391–407, 1955.

Hetényi, G., Cattin, R., Brunet, F., Bollinger, L., Vergne, J., Nábělek, J. L., and Diament, M.: Density distribution of the India plate beneath the Tibetan plateau: Geophysical and petrological constraints on the kinetics of lower-crustal eclogitization, Earth Planet. Sc. Lett., 264, 226–244, 2007.

Hilairet, N., Reynard, B., Wang, Y., Daniel, I., Merkel, S., Nishiyama, N., and Petitgirard, S.: High-pressure creep of ser-pentine, interseismic deformation, and initiation of subduction, Science, 318, 1910–1913, 2007.

Hirauchi, K.-I., Katayama, I., and Kouketsu, Y.: Semi-
brittle deformation of antigorite serpentinite under forearc
mantle wedge conditions, J. Struct. Geol., 140, 104151,
https://doi.org/10.1016/j.jsg.2020.104151, 2020.

Hirth, G. and Guillot, S.: Rheology and tectonic significance of ser-
pentinite, Elements, 9, 107–113, 2013.

Hirth, G. and Kohlstedt, D.: Rheology of the upper mantle and the mantle wedge: A view from the experimentalists, Geophysical Monograph-American Geophysical Union, 138, 83–106, 2003.

Holland, T. and Powell, R.: An internally consistent thermodynamic
data set for phases of petrological interest, J. Metamorph. Geol.,
16, 309–343, 1998.

Idrissi, H., Bollinger, C., Boioli, F., Schryvers, D., and Cordier, P.: Low-temperature plasticity of olivine revisited with in situ TEM nanomechanical testing, Sci. Adv., 2, e1501671, https://doi.org/10.1126/sciadv.1501671, 2016.

Idrissi, H., Samaee, V., Lumbeeck, G., van der Werf, T., Pardoen, T., Schryvers, D., and Cordier, P.: In situ quantitative tensile testing of antigorite in a transmission electron microscope, Journal of Geophysical Research: Solid Earth, 125, e2019JB018 383, 2020.

Jammes, S. and Huismans, R. S.: Structural styles of mountain building: Controls of lithospheric rheologic stratification and ex-tensional inheritance, J. Geophys. Res.-Sol. Ea., 117, B10403, https://doi.org/10.1029/2012JB009376, 2012.

Jammes, S., Manatschal, G., Lavier, L., and Masini, E.: Tectonosed-
imentary evolution related to extreme crustal thinning ahead of a
propagating ocean: Example of the western Pyrenees, Tectonics,
28, TC4012, https://doi.org/10.1029/2008TC002406, 2009.

Jammes, S., Huismans, R. S., and Muñoz, J. A.: Lateral variation in structural style of mountain building: controls of rheological and rift inheritance, Terra Nova, 26, 201–207, 2014.

Jaquet, Y. and Schmalholz, S. M.: Spontaneous ductile crustal shear zone formation by thermal softening and related stress, temper-ature and strain rate evolution, Tectonophysics, 746, 384–397, 2018.

Jaquet, Y., Duretz, T., Grujic, D., Masson, H., and Schmalholz, S. M.: Formation of orogenic wedges and crustal shear zones by thermal softening, associated topographic evolution and application to natural orogens, Tectonophysics, 746, 512–529, 2018.

Kameyama, M., Yuen, D. A., and Karato, S.-I.: Thermal-mechanical effects of low-temperature plasticity (the Peierls mechanism) on the deformation of a viscoelastic shear zone, Earth Planet. Sc. Lett., 168, 159–172, 1999.

Kaus, B. J., Connolly, J. A., Podladchikov, Y. Y., and Schmalholz, S. M.: Effect of mineral phase transitions on sedimentary basin subsidence and uplift, Earth Planet. Sc. Lett., 233, 213–228, 2005.

Kiss, D., Candioti, L. G., Duretz, T., and Schmalholz, S. M.: Thermal softening induced subduction initiation at a passive margin, Geophys. J. Int., 220, 2068–2073, 2020.

Kissling, E.: Deep structure of the Alps – what do we really know?, Phys. Earth Planet. In., 79, 87–112, 1993.

Kissling, E. and Schlunegger, F.: Rollback orogeny model for the evolution of the Swiss Alps, Tectonics, 37, 1097–1115, 2018.

Kronenberg, A. K., Kirby, S. H., and Pinkston, J.: Basal slip and mechanical anisotropy of biotite, J. Geophys. Res.-Sol. Ea., 95, 19257–19278, 1990.

Lamb, S. and Davis, P.: Cenozoic climate change as a possible cause for the rise of the Andes, Nature, 425, 792–797, 2003.

Lardeaux, J.-M.: Deciphering orogeny: a metamorphic perspec-
tive. Examples from European Alpine and Variscan belts: Part
I: Alpine metamorphism in the western Alps. A review, B. Soc.
Géol. Fr., 185, 93–114, 2014.

Le Breton, E., Brune, S., Ustaszewski, K., Zahirovic, S., Seton, M., and Müller, R. D.: Kinematics and extent of the Piemont-Liguria Basin – implications for subduction processes in the Alps, Solid Earth, 12, 885–913, https://doi.org/10.5194/se-12-885-2021, 2021.

Li, Z. and Gerya, T. V.: Polyphase formation and exhumation of high-to ultrahigh-pressure rocks in continental subduction zone: Numerical modeling and application to the Sulu ultrahigh-pressure terrane in eastern China, J. Geophys. Res.-Sol. Ea., 114, B09406, https://doi.org/10.1029/2008JB005935, 2009.

Mackwell, S., Zimmerman, M., and Kohlstedt, D.: High-
temperature deformation of dry diabase with application to tec-
tonics on Venus, J. Geophys. Res.-Sol. Ea., 103, 975–984, 1998.

Malavieille, J.: Impact of erosion, sedimentation, and structural her-
itage on the structure and kinematics of orogenic wedges: Analog
models and case studies, GSA Today, 20, 4–10, 2010.

Malinverno, A. and Ryan, W. B.: Extension in the Tyrrhenian Sea and shortening in the Apennines as result of arc migration driven by sinking of the lithosphere, Tectonics, 5, 227–245, 1986.

Malusà, M. G., Faccenna, C., Baldwin, S. L., Fitzgerald, P. G., Ros-
setti, F., Balestrieri, M. L., Danišík, M., Ellero, A., Ottria, G.,
and Piromallo, C.: Contrasting styles of (U) HP rock exhumation
along the Cenozoic Adria-Europe plate boundary (Western Alps,
Calabria, Corsica), Geochem. Geophy. Geosy., 16, 1786–1824, 2015.

Malvoisin, B., Austrheim, H., Hetényi, G., Reynes, J., Hermann, J., Baumgartner, L. P., and Podladchikov, Y. Y.: Sustainable densification of the deep crust, Geology, 48, 673–677, 2020.

Manatschal, G. and Müntener, O.: A type sequence across an ancient magma-poor ocean–continent transition: the example of the western Alpine Tethys ophiolites, Tectonophysics, 473, 4–19, 2009.

Mancktelow, N. S. and Pennacchioni, G.: Why calcite can be stronger than quartz, J. Geophys. Res.-Sol. Ea., 115, B01402, https://doi.org/10.1029/2009JB006526, 2010.

Manzotti, P., Bellevre, M., Zucali, M., Robyr, M., and Engi, M.: The tectonometamorphic evolution of the Sesia–Dent Blanche nappes (internal Western Alps): review and synthesis, Swiss J. Geosci., 107, 309–336, 2014.

McCarthy, A., Chelle-Michou, C., Müntener, O., Arculus, R., and Blundy, J.: Subduction initiation without magmatism: The case of the missing Alpine magmatic arc, Geology, 46, 1059–1062, 2018.

McCarthy, A., Tugend, J., Mohn, G., Candioti, L., Chelle-Michou, C., Arculus, R., Schmalholz, S. M., and Müntener, O.: A case of Ampferer-type subduction and consequences for the Alps and the Pyrenees, Am. J. Sci., 320, 313–372, 2020.

Mohn, G., Manatschal, G., Beltrando, M., and Haupert, I.: The role of rift-inherited hyper-extension in Alpine-type orogens, Terra Nova, 26, 347–353, 2014.

Muñoz, J. A.: Evolution of a continental collision belt: ECORS-Pyrenees crustal balanced cross-section, in: Thrust tectonics, 235–246, Springer, 1992.

Nicolas, A., Hirn, A., Nicolich, R., and Polino, R.: Lithospheric wedging in the western Alps inferred from the ECORS-CROP traverse, Geology, 18, 587–590, 1990.

Pelletier, L., Müntener, O., Kalt, A., Vennemann, T. W., and Belgya, T.: Emplacement of ultramafic rocks into the continental crust monitored by light and other trace elements: An example from the Geisspfad body (Swiss-Italian Alps), Chem. Geol., 255, 143–159, 2008.

Petri, B., Duretz, T., Mohn, G., Schmalholz, S. M., Karner, G. D., and Müntener, O.: Thinning mechanisms of heterogeneous continental lithosphere, Earth Planet. Sc. Lett., 512, 147–162, 2019.

Platt, J.: Dynamics of orogenic wedges and the uplift of high-pressure metamorphic rocks, Geol. Soc. Am. Bull., 97, 1037–1053, 1986.

Popov, A. and Sobolev, S.: SLIM3D: A tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology, Phys. Earth Planet. In., 171, 55–75, 2008.

Raimbourg, H., Jolivet, L., and Leroy, Y.: Consequences of progressive eclogitization on crustal exhumation, a mechanical study, Geophys. J. Int., 168, 379–401, 2007.

Raleigh, C. B. and Paterson, M.: Experimental deformation of serpentinite and its tectonic implications, J. Geophys. Res., 70, 3965–3985, 1965.

Ramberg, H.: Gravity, deformation and the earth’s crust: in theory, experiments and geological application, Academic press, 1981.

Ranalli, G.: Rheology of the Earth, Springer Science & Business Media, 1995.

Rubie, D. C.: The catalysis of mineral reactions by water and restrictions on the presence of aqueous fluid during metamorphism, Mineral. Mag., 50, 399–415, 1986.

Ruh, J. B., Kaus, B. J., and Burg, J.-P.: Numerical investigation of deformation mechanics in fold-and-thrust belts: Influence of rheology of single and multiple décollements, Tectonics, 31, TC3005, https://doi.org/10.1029/2011TC003047, 2012.

Rummel, L., Baumann, T. S., and Kaus, B. J.: An autonomous petrological database for geodynamic simulations of magmatic systems, Geophys. J. Int., 223, 1820–1836, 2020.

Rybacki, E. and Dresen, G.: Deformation mechanism maps for feldspar rocks, Tectonophysics, 382, 173–187, 2004.

Schenker, F. L., Schmalholz, S. M., Moulas, E., Pleuger, J., Baumgartner, L. P., Podladchikov, Y., Vrijmoed, J., Buchs, N., and Müntener, O.: Current challenges for explaining (ultra) high-pressure tectonism in the Pennine domain of the Central and Western Alps, J. Metamorph. Geol., 33, 869–886, 2015.

Schierjott, J., Rozel, A., and Tackley, P.: On the self-regulating effect of grain size evolution in mantle convection models: application to thermochemical piles, Solid Earth, 11, 959–982, https://doi.org/10.5194/se-11-959-2020, 2020.

Schmalholz, S., Podladchikov, Y., and Schmid, D.: A spectral/finite difference method for simulating large deformations of heterogeneous, viscoelastic materials, Geophys. J. Int., 145, 199–208, 2001.

Schmalholz, S. M. and Fletcher, R. C.: The exponential flow law applied to necking and folding of a ductile layer, Geophys. J. Int., 184, 83–89, 2011.

Schmalholz, S. M., Medvedev, S., Lechmann, S. M., and Podladchikov, Y.: Relationship between tectonic overpressure, deviatoric stress, driving force, isostasy and gravitational potential energy, Geophys. J. Int., 197, 680–696, 2014.

Schmalholz, S. M., Duretz, T., Hetényi, G., and Medvedev, S.: Distribution and magnitude of stress due to lateral variation of gravitational potential energy between Indian lowland and Tibetan plateau, Geophys. J. Int., 216, 1313–1333, 2019.

Schmid, S., Boland, J., and Paterson, M.: Superplastic flow in fine-grained limestone, Tectonophysics, 43, 257–291, 1977.

Schmid, S. M., Kissling, E., Diehl, T., van Hinsbergen, D. J., and Molli, G.: Ivrea mantle wedge, arc of the Western Alps, and kinematic evolution of the Alps–Apennines orogenic system, Swiss J. Geosci., 110, 581–612, 2017.

Shreve, R. L. and Cloos, M.: Dynamics of sediment subduction, melange formation, and prism accretion, J. Geophys. Res.-Sol. Ea., 91, 10229–10245, 1986.

Simpson, G. D.: Mechanical modelling of folding versus faulting in brittle–ductile wedges, J. Struct. Geol., 31, 369–381, 2009.

Sizova, E., Gerya, T., and Brown, M.: Contrasting styles of Phanerozoic and Precambrian continental collision, Gondwana Res., 25, 522–545, 2014.

Spitz, R., Bauville, A., Epard, J.-L., Kaus, B. J. P., Popov, A. A., and Schmalholz, S. M.: Control of 3-D tectonic inheritance on fold-and-thrust belts: insights from 3-D numerical models and application to the Helvetic nappe system, Solid Earth, 11, 999–1026, https://doi.org/10.5194/se-11-999-2020, 2020.

Stern, R. J.: Subduction initiation: spontaneous and induced, Earth Planet. Sc. Lett., 226, 275–292, 2004.

Stern, R. J. and Gerya, T.: Subduction initiation in nature and models: A review, Tectonophysics, 746, 173–198, 2018.
Stixrude, L. and Lithgow-Bertelloni, C.: Thermodynamics of mantle minerals-II. Phase equilibria, Geophys. J. Int., 184, 1180–1213, 2011.

Sutra, E., Manatschal, G., Mohn, G., and Unternehr, P.: Quantification and restoration of extensional deformation along the Western Iberia and Newfoundland rifted margins, Geochem. Geophy. Geosy., 14, 2575–2597, 2013.

Teixell, A., Labaume, P., Ayarza, P., Espurt, N., de Saint Blanquat, M., and Lagabrielle, Y.: Crustal structure and evolution of the Pyrenean-Cantabrian belt: A review and new interpretations from recent concepts and data, Tectonophysics, 724, 146–170, 2018.

Thielmann, M. and Kaus, B. J.: Shear heating induced lithospheric-scale localization: Does it result in subduction?, Earth Planet. Sc. Lett., 359, 1–13, 2012.

Toussaint, G., Burov, E., and Jolivet, L.: Continental plate collision: Unstable vs. stable slab dynamics, Geology, 32, 33–36, 2004.

Turcotte, D. and Schubert, G.: Geodynamics, Cambridge University Press, 2014.

van Hunen, J., van den Berg, A. P., and Vlaar, N. J.: Latent heat effects of the major mantle phase transitions on low-angle subduction, Earth Planet. Sc. Lett., 190, 125–135, 2001.

Vanderhaeghe, O., Medvedev, S., Fullsack, P., Beaumont, C., and Jamieson, R. A.: Evolution of orogenic wedges and continental plateaux: insights from crustal thermal–mechanical models overlying subducting mantle lithosphere, Geophys. J. Int., 153, 27–51, 2003.

Warren, C. J., Beaumont, C., and Jamieson, R. A.: Formation and exhumation of ultra-high-pressure rocks during continental collision: Role of detachment in the subduction channel, Geochem. Geophy. Geosy., 9, Q04019, https://doi.org/10.1029/2007GC001839, 2008.

Weijermars, R. and Schmeling, H.: Scaling of Newtonian and non-Newtonian fluid dynamics without inertia for quantitative modelling of rock flow due to gravity (including the concept of rheological similarity), Phys. Earth Planet. In., 43, 316–330, 1986.

Willett, S., Beaumont, C., and Fullsack, P.: Mechanical model for the tectonics of doubly vergent compressional orogens, Geology, 21, 371–374, 1993.

Willett, S. D.: Orogeny and orography: The effects of erosion on the structure of mountain belts, J. Geophys. Res.-Sol. Ea., 104, 28957–28981, 1999.

Wilson, J. T.: A new class of faults and their bearing on continental drift, Nature, 207, 343–347, 1965.

Wilson, R., Houseman, G., Buiter, S., McCaffrey, K., and Doré, A.: Fifty years of the Wilson Cycle concept in plate tectonics: an overview, Geol. Soc. Lond. Spec. Publ., 470, 1–17, 2019.

Winter, J. D.: Principles of igneous and metamorphic petrology, Pearson education, 2013.

Workman, R. K. and Hart, S. R.: Major and trace element composition of the depleted MORB mantle (DMM), Earth Planet. Sc. Lett., 231, 53–72, 2005.

Yamato, P., Agard, P., Burov, E., Le Pourhiet, L., Jolivet, L., and Tiberi, C.: Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural P-T-t data (Schistes Lustrés, western Alps), J. Geophys. Res.-Sol. Ea., 112, B07410, https://doi.org/10.1029/2006JB004441, 2007.

Yamato, P., Duretz, T., May, D. A., and Tartese, R.: Quantifying magma segregation in dykes, Tectonophysics, 660, 132–147, 2015.

Yamato, P., Duretz, T., and Angiboust, S.: Brittle/ductile deformation of eclogites: insights from numerical models, Geochem. Geophy. Geosy., 20, 3116–3133, https://doi.org/10.1029/2019GC008249, 2019.

Yang, J., Lu, G., Liu, T., Li, Y., Wang, K., Wang, X., Sun, B., Faccenda, M., and Zhao, L.: Amagmatic subduction produced by mantle serpentinization and oceanic crust delamination, Geophys. Res. Lett., 47, e2019GL086257, https://doi.org/10.1029/2019GL086257, 2020.

Zhao, L., Paul, A., Guillot, S., Solarino, S., Malusà, M. G., Zheng, T., Aubert, C., Salimbeni, S., Dumont, T., Schwartz, S., Zhu, R., and Wang, Q.: First seismic evidence for continental subduction beneath the Western Alps, Geology, 43, 815–818, 2015.

Zhao, L., Malusà, M. G., Yuan, H., Paul, A., Guillot, S., Lu, Y., Stehly, L., Solarino, S., Eva, E., Lu, G., Bodin, T., CIFALPS Group, and AlpArray Working Group: Evidence for a serpentinized plate interface favouring continental subduction, Nat. Commun., 11, 1–8, 2020.