Plume —lid interactions during the Archean and implications for the generation of early continental crust

Many Archean terranes are interpreted to have a tectonic and metamorphic evolution that indicates intra-crustal reorganization driven by lithospheric-scale gravitational instabilities. These processes are associated with the production of a signi cant amount of felsic and ma c crust, and are widely regarded to be a consequence of plume-lithosphere interactions. The juvenile Archean felsic crust is made predominantly of rocks of the tonalite–trondhjemite–granodiorite (TTG) suite, which are the result of partial melting of hydrous metabasalts. The geodynamic processes that have assisted the production of juvenile felsic crust, are still not well understood. Here, we perform 2D and 3D numerical simulations coupled with the state-of-the-art of petrological thermodynamical modelling to study the tectonic evolution of a primitive Archean oceanic plateau with particular regard on the condition of extraction of felsic melts. In our numerical simulations, the continuous emplacement of new, dry ma c intrusions and the extraction of the felsic melts, generate an unstable lower crust which drips into the mantle soon after the plume arrival. The subsequent tectonic evolution depends on the asthenosphere TP . If the TP is high enough (≥ 1500 ◦C) the entire oceanic crust is recycled within 2 Myrs. By contrast at low TP , the thin oceanic plateau slowly propagates generating plate-boundary like features.

A F 21 M , 2021 is the velocity vector component along the direction (i.e x,y,z). In order to incorporate the compaction and 60 expansion e ects due to the melt extraction, we introduced a source/sink (S, see Tab. S1 and below for further 61 explanations) term in the non-linear residual calculations. S has an indirect e ect on the momentum equation, because 62 it a ects the pressure and velocities. Therefore, the momentum equation retains its canonical form: 64 where is the deviatoric stress tensor, P is the pressure, is the density and is the gravity vector component 65 along the i direction. 66 LaMEM employs viscous-elasto-plastic rheological constitutive equations that connect the deviatoric stress tensor Wherė is the total strain rate tensor, the superscript el,vis,pl indicate the elastic (el), viscous (vis) and plastic 72 (pl) strain rate. ⋄ is the Jaumann objective stress rate, G is the shear modulus. is the spin tensor,̇ is the plastic 73 multiplier, Q is the plastic ow potential, which is equal to the second invariant of the deviatoric stress tensor ( ), 74 according to a dilatation-free non associative Drucker-Prager ow rule. is the viscosity which is computed using is the pre-exponential factor, n is the stress exponent,̇ is the second invariant of the strain rate tensor, R is the 78 gas constant, E and V are respectively the activation energy and volume.

79
The plastic rheology is enforced by the application of the Drucker-Prager failure criterion (Drucker & Prager, 1952): where is the heat capacity, / is the objective temperature time derivative, k is the heat conductivity, is the 93 radiogenic heat productivity ( = , where A is the amount of heat produced per unit of mass). H is the adiabatic 94 heating and is computed according the eq. 9 (in which is the thermal expansivity) while H is the dissipative heating. 95 We neglect latent heat of melting because its e ect is negligible compared with shear and adiabatic heating.
118 is the total amount of melt extracted from a rock type along z direction at xed x,y coordinates; is 119 the volumetric correction, which on average is 0.9; k, is the k ℎ node along z direction; is the last node of the 120 numerical domain along z direction.

121
The total amount of melt extracted is converted into extrusive and intrusive igneous rocks: 123 Where Vol and Vol is the extruded (e usive) and intruded volume, respectively. I represents the volumetric 124 amount of intrusion with respect to the total volume of melt extracted. The e usive volcanic rocks are treated as if 125 they were sediments: is converted into an e ective thickness dividing it for the current nite di erence cell 126 area, then is used to shift the (internal) free surface (see Fig. S1).

127
The intrusions are emplaced in a speci c area of the crust (see The source term is computed assuming that that extraction and injection locally shrink or expand the control volume 136 during the timestep. It is an e ective volumetric deformation (1/ ).  and M , and the long-term evolution of the system. 181 We mainly focus on the and/or conditions at which both felsic and ma c melts are extracted because this 182 information can be used as a compositional proxy. We classify the felsic melts extracted into two main categories:

183
TTG-like melts (TTG-opt) and non-TTG-like melts. TTG-like melts are the melts extracted within a -range that 184 produces a bulk composition that best ts the chemistry of the natural TTG (TTG-opt). We further split TTG-like 185 melts into three categories following Moyen (2011) as a function of the pressure of extraction within the optimum eld: 186 low-pressure TTGs (LP-TTG, 0.6-1.0 GPa, 800-1000 • C), medium-pressure TTGs (MP-TTG, 1.0-1.8 GPa, 800-1000 • C) 187 and high-pressure TTGs (HP-TTGs,>1.8 GPa, 800-1000 • C). Non-TTG-like melts are all other felsic melts of a broadly 188 granitic composition that did not form in the optimal P-T eld identi ed above. Ma c melts are tracked using their 189 extraction temperature, which can give insights on mantle source (plume or asthenosphere or contaminated mantle).

190
In order to describe the macroscopic structural features, we introduce the following terminology: external zone 191 represents the initial lithosphere (sometimes referred as old-lithosphere); internal zone is the newly produced 192 composite crust produced during the plume-lithosphere interaction and a ected by RTIs; boundary zone is the limit 193 between internal and external zone, dominated by the delamination of the lower crust of the external zone. (<1500 • C), we assumed that the temperature of the lower boundary is higher than the value expected for the given 233 to increases the convective vigors. The initial geothermal gradient of the lithosphere is a linear double stage geotherm, 234 featuring a T ℎ of 800 • C and a temperature at the base of the lithosphere that depends on the asthenospheric T .   3.3 E ect of the potential temperature of the asthenosphere 298 We explore the e ect of the T of the asthenosphere using 3D (run for 10 Myrs) and 2D (run for 40 Myrs) numerical 299 experiments. Our results can be divided into two main regimes: high mantle potential temperature (high T ) and low 300 mantle potential temperature regimes (low ).

354
The experiment featuring higher mantle potential temperature ( = 1550 • ) shows the same evolution (see Fig 6).

355
The only di erence that can be reported is related to the absolute quantity of new crust which is higher than the 356 reference scenarios. We additionally tested the e ect of changing the interval and the relative depth at which the 357 ma c intrusions are emplaced. The main geodynamic scenario is not changing but the pressure at which felsic melts 358 are extracted is higher than the reference scenarios resulting in larger amounts of TTG-like and MP-TTG-like melts.  The 2D experiment (see Fig 3) shows rst-order similarities with the 3D numerical experiment, but there are some  The experiment does not reach the quiescent stage (Fig 7, 8.524 Myrs). And the processes described above stabilize 386 the boundary zone generating a plate-like boundary feature. At the end of the time span explored, the crust can 387 be subdivided into three areas: the internal zone whose crust is mostly ma c with felsic intrusive bodies that are 388 emerging locally due to the RTIs (see Fig. 9, a, b d); an annular continental terrain at the rim of the internal zone 389 over-thrusting the external zone crust; a annular foreland basin whose depth is −1.0 km (Fig.9 c). The annular felsic 390 terrains have a topography that spans from 2.0 -3.0 km (Fig.9 c) and internal zone featuring a relative high topography 391 (0.6-1.5 km). The topographic gradient imposed by the annular felsic terrain is high and facilitates the over-thrusting 392 of the external zone ma c crust due to the gradient of gravitational potential energy between the annular felsic 393 terranes and fore-land basin (similar to the results of Rey et al. (2014)). In the internal zone the topographic heights 394 are associated with the amount of felsic crust, and the pattern of small basins is like those that have been observed in 395 the reference high scenarios (Fig.9 c-d).

396
The amount of new ma c crust generated increases linearly with time (see Fig 10,    associated with gravitational instabilities, that in case of a long-lasting plume are expected to be more e ective. 464 We introduce a novel approach to deal with the compositional variation due to magmatic processes; however, it  hybridization in our simulations, this thermal exchange must result in at least one of these processes (Bédard, 2006).

489
The mantle may be hydrated and contaminated by the dripped crustal material. These processes are feasible and    Our results show that several plate-like processes occur in a geodynamic scenario dominated by vertical displacements.

599
As noted by Sizova et al. (2015), the emergence of plate-like processes is most likely due to the compositional evolution 600 of the oceanic crust, which also has direct implications for the rock types diagnostic of subduction -such as blueschists

601
-that may form and be preserved in geological terranes through time  can evolve into modern-style subduction zones. If these zones connect and become a global network, the transition 606 from a stagnant lid regime to a mobile lid (plate tectonic) geodynamic regime may proceed.     . Column a) features sub-plots represents the compositional eld, while column b) shows several snapshots of the buoyancy contrast that controls the evolution of the numerical experiment, the black thick lines represent the internal free surface and the Moho topography, while the grey arrows represent the streamline of the velocity eld. Δ is computed by subtracting the density eld along x direction for a given depth from the average density of the considered depth (Δ ( , z) = ( , z) − (z), where x represents the x coordinate, and z represents the given depth). Therefore, this parameter represents the excess (negative values) and the de ciency (positive) of mass with respect to the background density. The snapshots of the 3rd stage represents the full numerical domain during the quiescent stage. c) represents the full compositional domain, while d) represents both the volumetric fraction of melt extracted from the mantle during the whole simulation, and the volumetric fraction of felsic crust components within the crust. The black thick lines represent the free surface and the Moho depth.  Fig. 4 are highlighted in red. The last row is the 2D numerical experiments, which feature the same input parameter of the high reference scenarios. The evolution of this experiment is shown in Fig. 3 .a) Felsic melt extraction conditions: left axis and black thick line represents the cumulative volume of felsic melts produced during the simulation (i.e. the felsic crustal growth curve); the right axis and the red thick line represents the cumulative volumetric fraction of TTG-like melts produced during the simulation. The shaded green, blue and red area represents LP-TTG, MP-TTG and HP-TTG like melts cumulative volumetric fraction respectively. b) Ma c melts extraction conditions: left axis and black thick line represents the cumulative volume of ma c melts extracted (i.e. the ma c crustal growth curve) during the simulation; right axis and red thick line represent the temperature of extraction of ma c melts ( ).The blue shaded area represents the temperature conditions in which 5-25% and 75-95% of melts are extracted. While the red shaded area represents the standard deviation of the during the simulation evolution. The initial temperature and nal temperature are highlighted in the right of each plot to highlight the e ective cooling of the mantle during the whole simulation. However, the e ective thickness is de ned in the cell center of the main numerical grid. Therefore this values must 661 be linearly interpolated into the main grid. After the interpolation, the free surface is shifted accordingly to the local 662 and the area between the previously computed free surface and the new one is lled with e usive rocks (see S14).

663
The cumulative melt production is computed by summing the total melt extracted along z direction, at every timestep:  Crust and lithosphere thickness , F n.d. Average amount of felsic crust components along z direction and proportion of felsic crust components at x,y,z coordinates.

/ [ 3 ]
Total amount of ma c or felsic crust produced. Table S2: List of the petrophysical properties and phase diagrams. The symbol and relative unit of measures are listed Tab. S1. All rock types share the same: shear modulus, ( = 40 ); initial and nal friction angle ( 0 = 30 • and 1 = 9 • ); initial and nal cohesion ( 0 = 10 and 1 = 3 MPa); heat capacity ( = 1200 / / ); thermal conductivity ( = 3). are the normalized total volume of low, middle and high pressure TTG-like melts. , and are the volumetric weighted proportion of felsic melts produced within low pressure conditions (≤ 1.0 GPa), middle pressure condition (1.0≥P≤1.8 GPa) and high pressure ( 1.8 GPa). , and represents the volumetric weighted proportion of felsic melts produced at low temperature condition (≤ 800 • C), middle temperature condition (800-1000 • )C and high temperature condition (≥ 1000 • C) respectively. and are the cumulative amount of felsic and ma c melts produced within the experiments. The resolution tests lists only the nal proportion of the TTG, and the nal cumulative volume of both ma c and felsic melts and the nal average crustal thickness ( ). The input parameters are the same of the Test3L, except for the number of nodes employed, that results into a di erent dx and dz and nite di erence cell area (A) A F 21 M , 2021  Figure S1: Simple representation of the free surface. 1) All the velocity's component are interpolated from the main grid node to the actual z coordinate of the free surface; 2) The z component is multiplied for the actual time-step (dt), then the free surface coordinate are shifted accordingly to this displacement; 3) The total volume of e usion is de ned in the cell center. To retrieve the e ective thickness is divided by the cell area;4) The free surface is again shifted according to and the di erence between the previous surface and the new one is lled with the e usive phase (e.g. BS0%).