Numerical modelling
The presented thermomechanical models were obtained by solving the conservation equation for a steady state momentum, transient heat conservation and incompressible mass conservation equations:
where v is the velocity vector, T is the temperature, k is the thermal conductivity, ρ is the density, cp is the heat capacity, Qr is the radiogenic heat production, τ is the deviatoric stress tensor, is the deviatoric strain rate tensor, P is the pressure and g is the gravity acceleration vector. The term describes the production of heat by viscoplastic dissipation (shear heating).
The density field evolves according the following equation of state:
whereρ0 is the reference density, α is the thermal expansivity, β is the compressibility, T0 and P0 are the reference temperature and pressure which were respectively set to 0 C and 105 Pa.
The effective viscosity () relates the deviator stress and strain rate tensor in the following fashion:
and is computed in order to satisfy a viscoelastoplastic rheological model:
where the v, e, and p superscripts correspond to viscous, elastic and plastic portions and the surperscripts dis and Peierls refers to the dislocation and Peierls creep mechanisms.
The viscous strain rate is computed as:
where A is a prefactor, Q is the activation energy, n is the stress exponent R is the universal gas constant and f is a correction factor 43. The subscripts II stand for the square root of the second tensor invariant. For rheological parameters used in the reference model, see Table 1. The elastic strain rate is written as:
where G is the shear modulus (set to 1010 Pa).
The plastic strain rate takes the form of:
where is the friction angle and is the cohesion (for friction angle and cohesion values of the reference model see Table 1). We do not apply any plastic strain softening.
In the mantle lithosphere, the Peierls mechanism is also activated and its strain rate is computed as:
where the effective strain rate is spelled as:
where the parameters s is the effective stress exponent (Tdependent), QPeierls is the activation energy (=540 J/mol), σPeierls is the Peierls stress (=8.5.109 Pa), EPeierls (=5.7.1011 s1), q (=2.0), andγ (=0.1) 44. Peierls creep stress is computed using a regularised formulation 45.
The temperature is kept constant at both the upper (0oC) and lower boundaries (1330oC) and the heat flow is set to 0oC across the right and left boundaries. A plate convergence rate of 3 cm/year is achieved by prescribing constant normal inflow velocities of Vin=1.5 cm/year along the upper 140 km of the two model sides, while mass conservation is satisfied by gradually increasing outflow below 140 km. The shear stress is set to zero along the left, right and lower boundaries. The upper boundary is a true free surface that dynamically evolves with time 17.
The initial temperature field is obtained by solving the steady state heat equation (neglecting shear heating) using reference thermal parameters (Table 1), excepted below the lithospheric mantle where the conductivity was set artificially high in order to produce a quasiadiabatic asthenosphere. The initial topography is set to 0 km.
The conservation equations are discretized using a finite difference/markerincell technique 46. The global linearized systems of equations are solved using a directiterative method 47. Nonlinear iterations are used at both local and global levels. At the local level, Newton iterations ensure exact partitioning of strain rates and correct evaluation of effective viscosity 48,49. At the global level, Picard iterations are employed to bestsatisfy mechanical equilibrium equations (to an absolute tolerance of 106 and within a maximum of 20 iterations).
Table 1
Thermal and rheological parameters used for different compositions in the reference model. The heat capacity (Cp) and the compressibility (β), were set to 1,050 J kg− 1 K− 1 and 10− 11 Pa− 1 for all compositions, respectively. Rheological parameters (preexponential factor (A), stress exponent (n), and creep activation energy (Q)) are set according to flow laws of mica 50, Westerly granite 24, mafic granulite 25, Maryland diabase 23, dry olivine 51, and serpentinite 52, from top to bottom in the table. Other material properties: ρ = density, k = thermal conductivity, Qr = radiogenic heat production, α = coefficient of thermal expansion, C = cohesion, ϕ = friction angle.
ρ (kg.m− 3)

k (W.m− 1.K− 1)

Qr (W.m− 3)

α (K− 1)

C (MPa)

ϕ

A (Pan.s− 1)

n

Q (J.mol− 1)

Sedimentary cover (mica)

2700

2.55

2.9e6

3.0e5

10

15

1.0e138

18

51.0e3

Continental upper crust (Westerly granite)

2750

2.8

1.65e6

3.0e4

10

30

3.1623e26

3.3

186.5e3

Continental lower crust (mafic granulite)

2900

2.8

1.65e6

3.0e4

10

30

8.8334e22

4.2

445.0e3

Oceanic crust (Maryland diabase)

2900

3.0

1.0e10

3.0e5

10

30

3.2e20

3.0

276.0e3

Lithospheric mantle (dry olivine)

3300

3.0

1.0e10

3.0e5

10

30

1.1e16

3.5

530.0e3

Asthenosphere (dry olivine)

3300

3.0

1.0e10

3.0e5

10

30

1.1e16

3.5

530.0e3

Weak zone (serpentinite)

2900

3.0

1.0e10

3.0e5

0

30

4.4738e38

3.8

8.9e3

Model geometry
The computational domain is a cross section of 1330 × 410 km. The model resolution is 1 km in both directions. The initial compositional geometry is inspired by reconstructions of preobduction geodynamic settings. It contains an oceanic domain (660 km wide) with a spreading ridge and a tilted weak zone in the center that ensures leftdipping intraoceanic subduction initiation. The thermal structure of the oceanic lithosphere is calculated by applying a halfspace cooling age model from 1.5 Myr at the center to 50 Myr at the edges of the ocean. The oceanic crust is 6 km thick and is overlain by a layer of uppermost sedimentary cover that linearly thickens from the ridge (0 km) towards the right edge of the continental domain (3 km). The transition from the continent to the ocean is defined by a passive margin geometry where the continental upper and lower crust linearly thin from 30 km to 5 km over the distance of 200 km. The uppermost sedimentary cover layer has a constant 3 km thickness over the continental domain.
Data availability
Model output files that support the results of this study are stored in the data repository of Utrecht University and are available upon request.