Article
Methods
Numerical modelling method
The regional two-dimensional magmatic-thermomechanical model
simulates the subduction of an oceanic plate under another oceanic
plate followed by collision of incoming passive continental margin
with the oceanic arc. The governing equations of conservation of
momentum, mass and energy are solved with the use of the code
I2VIS^40 , based on conservative finite differences and a non-diffusive
marker-in-cell technique applied on a staggered non-uniform Eulerian
grid. The 4,000 km × 1,000 km numerical model domain is resolved
with a non-uniform 2,041 × 381 rectangular grid that provides a high-
est grid resolution of 1 km in the 1,500-km-wide and 200-km-thick
subduction area of the model. More than 70 million Lagrangian mark-
ers are used to carry material properties, temperature, water content
and melt-depletion, which are initially distributed over the very dense
randomly perturbed Lagrangian grid. These Lagrangian markers move
according to the computed velocity field. Visco-plastic non-Newtonian
rheologies are used in all experiments. The model includes creeping
flow with thermal and chemical buoyancy forces, accounts for eclogite,
olivine-spinel and spinel-perovskite phase transitions, and considers
the effects of adiabatic, shear, latent and radioactive heating. Complete
details of this method, allowing for its reproduction, are provided
elsewhere^40. It should also be noted that the performed high-resolution
numerical experiments are central processing unit (CPU)-time and
storage intensive and require long running time (up to several months)
owing to the small time-step size caused by the low viscosity of the
hot mantle.
Model geometry, initial and boundary conditions
We conducted a series of numerical experiments with different mantle
potential temperature (Tp) of 1,450 °C, 1,500 °C and 1,550 °C (that is,
ΔT = 150 °C, 200 °C and 250 °C hotter than the present-day values,
representative of Precambrian conditions^8 ,^10 ,^41 ,^42 ) and subducting plate
velocities of 5–10 cm yr−1 (Extended Data Table 1, Extended Data Fig. 1).
The 500-km-wide subducting oceanic basin is prescribed between
the left-hand side of the continent and the incipient subduction zone
(Extended Data Fig. 1). It is widely accepted that the Archaean oceanic
crust was much thicker than the modern one because of the extensive
melting of the relatively hot mantle (for example, see refs.^10 ,^43 ,^44 ). We
assumed that the oceanic crust is 25 km thick and consists of layer of
hydrothermally altered basalts (7 km thick) underlain by an 18-km-thick
layer of gabbroic rocks with rheology of wet quartzite and plagioclase^45 ,
respectively (Extended Data Table 1). In an additional series of experi-
ments, the thickness of the oceanic crust increases linearly from 20 km
to 30 km with mantle potential temperature increasing from 1,450 °C
to 1,550 °C (Extended Data Table 1).
The model continental crust is 40 km thick and composed of three
layers^46. The upper and middle crust (each 10 km thick) are felsic, with
rheology of wet quartzite, whereas the lower crust (20 km thick) is
mafic, with rheology of plagioclase^45. In an additional series of experi-
ments, the upper and middle crusts are 15 km thick and the lower crust
is 10 km thick (Extended Data Table 1).
The mantle is represented by anhydrous peridotite, which is initially
subjected to depth-dependent melt-depletion in accordance with the
mantle potential temperature using the melting model of Katz et al.^47
(see the MatLab code provided with this paper; see ‘Code availabil-
ity’ section). The initiation of subduction is enabled by prescribing
an initial weak zone (inclination angle 30°) with the rheology of wet
olivine^45 and low brittle/plastic strength (Extended Data Table 2). A
prescribed constant plate velocity of 5–10 cm yr−1 is defined within the
distant intraplate region of the left-hand side of the model and drives
the spontaneously bending oceanic slab to the right. To provide bet-
ter comparison with previous Archaean subduction models driven by
slab pull^24 ,^25 we also performed numerical experiments in which forced
convergence is switched off at 5–8 Myr after the beginning of subduc-
tion (Extended Data Fig. 8).
The initial thermal structure of the oceanic lithosphere is defined
by the plate cooling model with a plate thickness of 95 km (ref.^48 ). The
initial cooling age of the oceanic lithosphere is taken to be 40 Myr in
all experiments (Extended Data Table 1), which is similar to a previous
numerical modelling study that investigated changes in the style of
oceanic–continental^23 and intraoceanic subduction^22. In the conti-
nental crust, temperature changes linearly from 0 °C at the surface to
600 °C, at the Moho boundary. This Moho temperature is thus about
100–150 °C higher than that of the modern shields and platforms (for
example, refs.^49 ,^50 ). In the continental mantle lithosphere, tempera-
ture changes linearly from the Moho to 188 km depth for ΔT = 150 °C,
196 km for ΔT = 200 °C and 205 km for ΔT = 250 °C, where it reaches
the asthenospheric temperature subjected to adiabatic gradient of
0.5 K km−1. Comparison of the piecewise linear geotherm we used with
a steady-state continental geotherm that takes into account crustal
radioactive heating reveals their close correspondence. Results of
numerical experiments based on these two different geotherms also
do not show any substantial differences.
All mechanical boundary conditions are free-slip except the upper
and lower boundaries, which remain open (see Baitsch-Ghirardello
et al.^51 ). The top surface of the lithosphere is treated as an inter-
nal free surface^52 by using an overlying 22-km-thick layer with low
viscosity (10^18 Pa∙s) and density of 1 kg∙m−3 for air (<18 km) and
1,000 kg∙m−3 for seawater (>18 km). This upper boundary evolves
by erosion and sedimentation according to the following Eulerian
transport equation^23 ,^51 ,^53 :z
tvvz
tvv∂
∂=−∂
∂es zxes−−xe,where x and z are horizontal and vertical coordinates, respectively, zes
is the vertical position of the surface as a function of the horizontal
distance x, vz and vx are the vertical and horizontal components of
the material velocity vector at the surface, and vs and ve are the sedi-
mentation and erosion rates, respectively. The sedimentation and
erosion rates correspond to the following relations^51 : vs = 0 mm yr−1,
ve = ve0, when z < zsea, and vs = vs0; and ve = 0 mm yr−1, when z > zsea, where
ve0 = 0.3 mm yr−1 and vs0 = 0.03 mm yr−1 are the imposed erosion and
sedimentation rates, respectively, and zsea = 18 km is the sealevel pre-
scribed in the model. The surface slope, φmax, for the accumulated
sedimentary prism varies up to 35° (tanφmax = 0.7). In some models,
by regulating fluxes through the open upper and lower boundaries,
we prescribed gradual lowering of the reference level for the free ero-
sion–sedimentation surface with time to guarantee high quality for
the free surface approximation^52 at the arc-continent collision stage.
This caused gradual increase in the thickness of the water layer in
respective models (see Fig. 1a–d). Testing of this condition showed
that it does not much affect results of numerical experiments in term
of deep processes.Petrological model
The density of rocks varies with pressure (P) and temperature (T)
according to the equation:ρρPT,0=[1−αT(−Tβ 00 )][1+(PP−)],where ρ 0 is the standard density at P 0 = 1 MPa and T 0 = 298 K, and α
and β are the coefficients of thermal expansion and compressibility,
respectively (Extended Data Table 2).
The model takes into account transformations of olivine into wad-
sleite and ringwoodite^54 and into bridgmanite in the mantle^55. Eclogi-
tization of subducted basaltic and gabbroic crust is taken into account
by linear increasing the density of the crust with pressure from 0% to