Numerical modelling
The presented thermo-mechanical 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 visco-plastic 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 visco-elasto-plastic 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 pre-factor, 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 (T-dependent), QPeierls is the activation energy (=540 J/mol), σPeierls is the Peierls stress (=8.5.109 Pa), EPeierls (=5.7.1011 s-1), 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 quasi-adiabatic asthenosphere. The initial topography is set to 0 km.
The conservation equations are discretized using a finite difference/marker-in-cell technique 46. The global linearized systems of equations are solved using a direct-iterative method 47. Non-linear 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 best-satisfy mechanical equilibrium equations (to an absolute tolerance of 10-6 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 (pre-exponential 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 (Pa-n.s− 1)
|
n
|
Q (J.mol− 1)
|
Sedimentary cover (mica)
|
2700
|
2.55
|
2.9e-6
|
3.0e-5
|
10
|
15
|
1.0e-138
|
18
|
51.0e3
|
Continental upper crust (Westerly granite)
|
2750
|
2.8
|
1.65e-6
|
3.0e-4
|
10
|
30
|
3.1623e-26
|
3.3
|
186.5e3
|
Continental lower crust (mafic granulite)
|
2900
|
2.8
|
1.65e-6
|
3.0e-4
|
10
|
30
|
8.8334e-22
|
4.2
|
445.0e3
|
Oceanic crust (Maryland diabase)
|
2900
|
3.0
|
1.0e-10
|
3.0e-5
|
10
|
30
|
3.2e-20
|
3.0
|
276.0e3
|
Lithospheric mantle (dry olivine)
|
3300
|
3.0
|
1.0e-10
|
3.0e-5
|
10
|
30
|
1.1e-16
|
3.5
|
530.0e3
|
Asthenosphere (dry olivine)
|
3300
|
3.0
|
1.0e-10
|
3.0e-5
|
10
|
30
|
1.1e-16
|
3.5
|
530.0e3
|
Weak zone (serpentinite)
|
2900
|
3.0
|
1.0e-10
|
3.0e-5
|
0
|
30
|
4.4738e-38
|
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 pre-obduction geodynamic settings. It contains an oceanic domain (660 km wide) with a spreading ridge and a tilted weak zone in the center that ensures left-dipping intra-oceanic subduction initiation. The thermal structure of the oceanic lithosphere is calculated by applying a half-space 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.