Buoyancy versus shear forces in building orogenic wedges

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, FB/FD=ArF, which controls the mode of orogenic wedge formation: ArF≈0.5 causes thrust-sheet-dominated wedges, ArF≈0.75 causes minor wedge formation due to relamination of subducted crust below the upper plate, and ArF≈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 ArF⪅0.5 due to the absence of high-grade metamorphic rocks, whereas for the Alps ArF≈1 during exhumation of high-grade rocks and ArF⪅0.5 during the post-collisional stage.
In the models, FD increases during wedge growth and subduction and eventually reaches magnitudes (≈18 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.


2 Numerical model 2.1 Mathematical model and numerical algorithm As common 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 (no inertia) flow under gravity and the influence of boundary tractions. Brittle-plastic deformation is controlled 95 by a Drucker-Prager yield function. We do not apply any strain softening, such as frictional or viscous strain softening. We also do not consider grain size evolution. However, thermal softening is active, resulting from the conservation of energy and temperature-dependent effective viscosities. 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 deformation . Appendix A provides a detailed description 100 of the algorithm.

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 1x1 km in horizontal x and in vertical z direction, respectively. A stressfree surface (Duretz et al., 2016a) is set initially at z = 0 km and to allow for dynamic build-up of topography, the topmost 105 +20 km of the model domain are left free of material. We impose constant material inflow/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°C and 1612°C, respectively.
Viscous deformation is modelled leveraging a combination of rheological flow laws such as dislocation, diffusion and Peierls creep. The effective density of the materials is either directly calculated via a linearized equation of state or pre-computed based 110 on equilibrium phase-diagram sections for specific bulk-rock composition (see Appendix A). For more detailed justifications of model assumptions and explanations for the implementation of surface processes and 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. These elliptical inclusions represent structural and petrological units inherited from previous orogenic cycles. Weak and strong inclusions may represent 115 metasedimentary units and mafic intrusions, respectively. Below the upper crust an 8 km thick lower crust is employed. The initial Moho is set to z = −33 km and the mantle lithosphere extends down to z = −120 km. In total, 12 weak elliptical inclusions over two vertical levels are included into the mantle lithosphere. These inclusions represent weaknesses caused by, for example, more hydrated peridotites (see also Barnhoorn et al. (2010); Petri et al. (2019) for more details). The elliptical inclusions help generating more realistic and asymmetric margin geometries (Duretz et al., 2016b;Petri et al., 2019). We 120 include the upper mantle and transition zone down to z = −660 km. The distribution and emplacement of the inclusions is described in detail in appendix A. The difference between mantle lithosphere and upper mantle and transition zone is due to temperature only, i.e. all material parameters are the same.  Fperh Details on the solution models can be found in the solution_model.dat data file in Perple_X.

155
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, as well as the subduction dynamics are documented in the video supplement of models REF (Candioti, 2020a) and GC1 (Candioti,160 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 necking zone of the left margin is ≈ 20 km wide. In total, the width of this margin is ca. 80-90 km. The right margin is in total ca. 140-160 km wide including an ≈ 60-80 km wide necking zone. At the end of the cooling period (109 Myr, Fig. 4b), an ≈ 360 km wide basin has opened, which is floored by exhumed mantle material. Convection in the upper mantle stabilises the thermal and mechanical thickness of the lithosphere to ca. 90-100 km (no velocity glyphs in Fig. 4b, see also discussion in Candioti et al. (2020) for more detail). At 112 Myr subduction is initiated at the transition from the proximal to the distal 170 continental margin in the right model side (see Sutra et al. (2013)  upper crustal block of the overriding plate. In AS18, the serpentinite channel is intersected by the two colliding plates at z≈-30 km (Fig. 6e). In contrast to AS1, the continental upper crust is partly sheared off the subducting plate at z≈-25 km in AS18 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 200 density changes compared to LEOS-models (compare Fig. 3g to f). Similar to AS1, a serpentinite channel forms after subduction initiation in AC1 (Fig. 7a). 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. 7b & c). In contrast to AS1, exhumation of the subducted upper continental crust has not occurred until 163 Myr (compare Fig. 6c to Fig. 7c) in model history. At 170 Myr, the subducted continental crust in AC1 has been exhumed to the surface through the weak subduction interface (Fig. 7d). In contrast to AS1, 205 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. 7d). Figure 8a & b shows the density field of models AS1 and AC1, respectively, computed by the algorithm at 170 Myr in model history. Compared to AS1, in AC1 the density of the material in the wedge is much more variable. The density of the subducted crustal units in AC1 is up to ca. 100 kg m −3 higher than in AS1. The resulting buoyancy contrasts are much smaller leading to significantly less uplift of topography in AC1 compared to AS1. In AS1, the 210 maximum elevation of topography exceeds 10 km; in AC1 the maximum elevation is locally ≈ 7 km, but high topographies are in average ≈ 5 km. The width of the wedge is ≈350 km and ≈275 km in AS1 and AC1, respectively. Compared to AS18, only minor volume of serpentinite material has been sheared off the subducting slab in AC18 (compare Fig. 6e to

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 220 (Fig. 9a). Continental upper crust of the subducting plate is buried to z≈-120 km ( Fig. 9a & b). At ca. 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 surface along the weak subduction interface breaking through the continental upper crust of the overriding plate (Fig. 9c). The serpentinite material is surrounding the exhumed crustal material (Fig. 9d) 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 225 second slab detachment occurs at z≈-250 km (see Candioti (2020b)). In contrast to AC18, less crustal volume is subducted in GC18 (Fig. 9e). Instead, the continental upper crust is largely sheared-off the subducting lithosphere at z≈-40 km (Fig. 9f) and several thrust sheets form (Fig. 9f). Compared to GC1, the relatively stronger serpentinite in GC18 is largely subducted and sediments originally deposited in the trench are incorporated into a growing orogenic thrust wedge ( Fig. 9f-h).

Evolution of buoyancy and shear forces
In order to investigate the relative impact of buoyancy and shear forces on collision, we quantify the temporal evolution of the horizontal driving force per unit length (Fig. 10b

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. 11a). We refer here to these variations as different modes of orogenic wedge formation. These varying modes 260 depend mainly on the shear resistance of upper crust and serpentinite, and its buoyancy contrast to the surrounding mantle material. The magnitude of Ar F is used to classify the different modes of orogenic wedge formation.
For low shear resistance of serpentinites, Ar F ≈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. 11a) or (ii) a decreased buoyancy contrast (CD-265 models) and low shear resistance of the upper crust (GC1, see Fig. 11c).
(2) A channel-flow mode of exhumation occurs for a decreased buoyancy contrast (complex-EOS) and increased shear resistance of the upper crust (AC1, see Fig. 11b).
In models of high serpentinite shear resistance, also two different modes are observed: (3)

Buoyancy vs. shear forces controlling modes of orogenic wedge formation
The shear resistance of (i) serpentinites and (ii) the upper continental crust directly impacts on the shear forces (compare GC1 275 to AC1 in Fig. 10b). Determining the shear resistance, or effective viscosity, of the lithosphere deforming under geological time scales remains challenging (Burov et al., 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 280 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 (PT) conditions remains elusive (David et al., 2018;Hirauchi et al., 2020, & reference therein). 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 285 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 these mod-290 els. Instead, strong serpentinite material leads to detaching and thrusting of crustal material escaping subduction already at shallow depths. In natural settings, the effective strength of serpentinites may vary along the upper regions of the subduction plate in direction 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 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 precomputed density tables are based on calculated equilibrium phase diagrams. In our calculations we have assumed H 2 O saturation and therefore our system is open with respect to its H 2 O content, but closed with respect to other elements. This is a valid approximation for the km-scale that we are considering. However, the preservation of high-grade metamorphic rocks at the near-surface 305 environment indicates that the crustal metamorphic rocks are not always at thermodynamic equilibrium as it 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- 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 315 changes the mode of orogenic wedge formation modelled here (Fig. 11). Hence, it is important to (1) better constrain the mechanical strength of crustal rocks and serpentinites under natural deformation conditions and (2) to proceed towards more realistic density structures including metamorphic reactions in numerical models of collisional orogen formation.

Subduction initiation
We follow here the modelling approach by Candioti et al. (2020). The margin geometry and thermal structure, prior to con-320 vergence, is generated during a modelled rifting and cooling period. This way, the generated passive margin and marine basin system is modelled internally-consistent with respect to the subsequent convergence. Subduction initiation is horizontally forced (Stern, 2004;Stern and Gerya, 2018;Crameri et al., 2020) and a major lithosphere shear zone forms around the transition from the distal to the proximal margin (see also discussion Candioti et al. (2020)). The ad hoc parameterized layer of serpentinite is not relevant for subduction initiation, as subduction is initiated also in the reference model without serpentinite 325 (see Fig.4f). Instead, geometrical focusing of stresses below the margin together with thermal softening and a temperaturedependent 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 TN m −1 are necessary to initiate subduction for quartz-and feldspar-dominated upper crust, respec-subduction initiation at an idealised, 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 parameterisation 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 335 necessary for subduction initiation. However, the minimum value of F D required for subduction initiation in our models should be determined in future studies.

Model 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., 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 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)  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 criticised 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 . (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 370 from seismic tomography showing no significant exhumation of lower crust (Schmid et al., 2017). In our models, significant synconvergent exhumation of upper crust can occur by either diapirism or channel-flow and is enabled by spatially localized upper plate extension (Fig. 11) and not by significant erosion. Exhumation of lower crust does not occur in our model, which 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 375 exhumation mechanisms of (ultra)high-pressure rocks and (iii) the pressure and temperatures paths and associated exhumation velocities in a subsequent study and we will not further discuss these issues 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. 10d). In fact, in our models plate driving forces eventually reach again the magnitude necessary

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). The evolution of the Alpine orogeny, however, involves the entire lithosphere in wedging (Nicolas et al., 1990). Some studies extended the "classical" wedge model to the lithospheric scale 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, which is likely important when 400 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 shear resistance and density structure of material.
Our models, therefore, presumably provide a more realistic insight into orogen dynamics compared to orogenic wedge models considering only crustal deformation.
As mentioned in Section 4.1, the density structure impacts on the force balance and reduces the mean topographic elevation.

405
This is important, because orogenic wedge models applied to the Western Alps have been criticised 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 would reduce significantly the convergence velocities at the mature stage of orogenic Absolute values for shear and buoyancy forces reported here strongly depend on the amount and strength of crustal material 415 carried into the subduction zone. How much crustal material has been involved in Alpine subduction depends on the preorogenic 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 numerical resolution (1 km) and probably strongly impacts on the size of the modelled wedge. At higher numerical resolution, the size of the mechanical heterogeneities could be reduced and the modelled rifted 420 margins might be thinner, leading to a more realistic pre-orogenic crustal thickness generated 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 evolution, i.e. the relative increase of buoyancy forces with continuous subduction, would be most likely unchanged. Therefore, our models likely provide a representative model for the relative evolution of shear and buoyancy forces building orogenic wedges.

Conclusions
Our models show that upward-directed buoyancy forces, caused by subduction of continental crust, can be as high as the shear forces induced by far-field plate convergence and should, therefore, be considered in models of continent collision and associ-ated orogenic wedge formation. Parameters controlling the force balance and the mode of deformation and exhumation within the evolving orogenic wedge are: (1) The shear resistance of serpentinites, controlling the shear resistance of the subduction 430 interface, and of the continental upper crust, controlling the maximal depth of crustal subduction.
(2) The density structure of the subducted material. The density structure significantly impacts on the mean topographic elevation and leads to more realistic topography evolution.
The strength of serpentinites and of the upper crust, as well as the density structure exert a first-order control on the formation of orogenic wedges. Varying these three parameters only in our models generates a wide spectrum of different orogenic 435 wedge modes, including thrust-sheet dominated wedges, buoyancy flow dominated wedges, and minor wedge formation due to significant relamination of subducted crust below the upper plate. The spatial variation of these three parameters may explain many differences in deformation style observed in natural collisional orogens.
The increase of upward-directed buoyancy forces during orogenic wedge growth causes an increase of the required horizontal driving forces, which may cause a slow down, or "choking" of the associated subduction. This increase of horizontal driving 440 forces may cause horizontally-forced subduction initiation in other regions. Therefore, quantifying the magnitudes of 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.
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 (Candioti, 2020a) and the subduction 445 and collision stage of the numerical simulation GC1 (Candioti, 2020b).

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 where i and j are spatial indices and repeated inidices 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 where P (negative mean stress) is pressure, T is temperature, ρ 0 is the material density at reference temperature (T 0 ) and pressure (P 0 ), α and β are material parameters accounting for density changes due to thermal expansion and isothermal compression, respectively, ∆T = T−T 0 and ∆P = P−P 0 . 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 equilibria models is stored in look-up tables and density values are read in by the applied algorithm 460 during simulation runtime 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 where δ ij is the Kronecker-Delta, η eff is the effective viscosity,ε eff ij are the effective deviatoric strain rate tensor components, where G is shear modulus, ∆t is the time step, τ o ij are 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 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: 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, f H2O and r are material parameters. The diffusion creep viscosity is calculated as 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 regularised form of Kameyama et al. (1999) as where s is an effective stress exponent:

485
where A P , γ and σ P a flow law parameters. Brittle-plastic material failure is controlled by the Drucker-Prager yield function that 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 ij . In Eq. A4, the effective viscosity is computed either the quasi-harmonic average of the visco-elastic contributions as ) or is equal to the viscosity η pla calculated according to Eq. A14. Rigid body rotation is computed analytically: where denotes matrix transposition, R is the rotation matrix, θ is the rotation angle and ω ij denote vorticity tensor compo-500 nents. Heat transfer is included into the model by expressing the energy balance equation w.r.t. temperature as where c P is the specific heat capacity, D/Dt is the material derivative, k is thermal conductivity, H A = Tαv z gρ includes contributions from adiabatic processes assuming lithostatic pressure conditions, H D = τ 2 ij / 2η 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 a = 30 km and the semi minor axis 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 1st elliptical inclusion's center x C n at each vertical level is calculated as 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 C n+1 is then calculated as where R dx is a random spacing ensuring that the ellipses do not overlap each other. The z-coordinate of the elliptical inclusion's center z C m at the vertical level m is calculated as where z Moho = 33 km is the initial depth of the Moho. Finally, all particles 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 where x M and z M are the horizontal and vertical coordinate of the marker, respectively. All random numbers used here are 525 seeded at a value of 197 using the C-function srand. Choosing the above mentioned values yields an increased number of weak elliptical inclusions in the center of the domain. This yields localisation of deformation without an additional perturbation of the marker field in the center of the domain.

Appendix C: Buoyancy and driving forces
In this study we use F D = 2 ×τ avg II (x) as a representative for the driving force acting on the deforming system, whereτ avg II (x) is 530 calculated averaging the vertical integral of the second invariant of the deviatoric stress tensor The reader is referred to Candioti et al. (2020) for further detail.
For calculation of the buoyancy force per unit length 535 the density difference, ∆ρ, between the subducted material and the surrounding mantle is integrated over the area Ω of material subducted below z=-40 km. To obtain a significant value for the negative (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 subtracte from F B . As a measure for buoyancy or shear forces dominating orogen dynamics, we here 540 define Ar F = F B /F D which is comparable to the Argand number (England and McKenzie, 1982). Table A1. Physical parameters used in the numerical simulations M1-6.